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SUMMARY 


The  main  objective  of  the  Combustor  Design  Model  Evaluation  program  is  to  evaluate,  develop 
and  validate  quantitatively  accurate  advanced  analytical  methods  and  models  to  meet  the 
challenges  of  designing  advanced  combustors.  The  main  challenge  is  to  design  combustors,  that 
meet  stringent  requirements  of  combustor  aerothermal  performance,  emissions,  operability, 
exhaust  signature  and  durability,  in  a  cost-  and  time-effective  manner  through  greater  reliance  on 
analytical  design  methods. 

The  specific  objectives  of  the  program  are  to: 

•  Develop  and  validate  more  accurate  computer-based  analytical  models  for  combustor  flows. 

•  Design  benchmark  quality  experiments  and  use  the  data  from  advanced  diagnostic 
measurements  for  such  development  and  validation. 

•  Incorporate  the  analytical  models  into  a  next-generation  gas  turbine  combustor  design  system 
and  assess  its  capabilities. 

The  joint  velocity-scalar  probability  density  function  (pdf)  method  (or  joint  pdf  method  for  short) 
is  the  focus  of  the  model  development  and  validation  efforts  at  Allison  under  this  program.  The 
pdf  method  was  selected  as  the  most  promising  candidate  for  further  development  as  a 
combustor  design  tool  after  a  review  of  several  current  and  evolving  approaches.  The  pdf 
method  overcomes  the  main  deficiencies  of  the  current  models;  turbulent  transport  and 
turbulence/chemistry  interactions  appear  in  closed  form  in  the  pdf  method  and  need  not  be 
modeled.  Significant  previous  development  work  with  the  method  had  demonstrated  its 
advantages  over  conventional  models  and  its  feasibility  for  development  as  a  combustor  design 
tool.  However,  considerable  work  needed  to  be  done  in  order  for  the  method  to  be  used  for  the 
complex  swirling  recirculating  flow  fields  in  practical  combustors. 

Benchmark  experimental  configurations  were  designed  in  conjunction  with  Wright  Laboratory, 
Combustion  Branch  (WL/POSC)  personnel  and  University  of  Dayton  Research  Institute  (UDRI) 
personnel.  Measurements  were  made  by  UDRI  personnel  at  WL/POSC  under  a  separate 
program.  These  data  as  well  as  data  available  from  other  sources  (such  as  data  from  NASA  Hot 
Section  Technology  [HOST]  programs)  were  used  for  model  validation  in  the  current  program. 

Substantial  progress  has  been  made  under  the  current  program  in  developing  and  validating  the 
pdf  method  for  complex  flows  such  as  in  gas  turbine  combustors.  The  primary  areas  of 
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progress  have  been  in  the  development  of  the  models  for  the  turbulent  frequency  and  the 
development  of  the  algorithms  for  determining  the  mean  pressure  field  within  the  pdf 
calculations.  Other  areas  have  included  the  treatment  of  swirling  flows  and  treatment  near  solid 
walls.  These  developments  are  key  to  the  application  of  the  pdf  method  to  complex  flows  of 
practical  interest.  The  pdf  method  has  been  applied  to  variety  of  swirling,  recirculating, 
nonreacting  and  reacting  flows  under  this  program.  Results  from  pdf  calculations  are  in 
excellent  agreement  with  detailed  benchmark  experimental  data  for  temperature  and  turbulent 
velocity  correlations  up  to  fourth  order.  The  results  from  the  joint  pdf  method  have  been 
assessed  against  the  results  from  other  conventional  methods  as  well  as  the  scalar  pdf  method  (a 
subset  of  the  joint  pdf  method).  The  pdf  methods  have  been  shown  to  be  more  accurate  than 
conventional  turbulent  combustion  models.  A  baseline  design  system  based  on  the  pdf  methods 
has  been  developed.  Areas  of  further  development  to  fully  utilize  the  design  system  have  been 
identified  and  discussed.  All  the  program  objectives  have  been  achieved. 
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1.  INTRODUCTION 


Future  high  performance  gas  turbine  engines  are  placing  exacting  requirements  on  combustor 
design  in  regard  to  aerothermal  performance,  emissions,  operability,  exhaust  signature,  durability 
and  combustor  envelope.  To  meet  these  requirements  in  a  cost-  and  time-effective  manner, 
improved  combustor  design  systems  are  needed  that  are  based  on  fundamentally  sound  and 
quantitatively  accurate  physicochemical  models  of  combustion  processes.  The  accuracy  in 
modeling  these  processes  is  key  to  the  accurate  prediction  of  critical  combustion  performance 
characteristics  such  as  efficiency,  emissions,  exit  and  wall  temperatures,  lean  blowout,  stability, 
flashback,  autoignition,  etc.,  in  advanced  high  performance  and  low  emissions  combustor 
concepts. 

Current  combustor  design  systems  lack  the  quantitative  accuracy  needed  for  the  reliable 
analytical  predictions  of  key  combustor  operating  characteristics.  This  deficiency  results  in  high 
costs  and  time  involved  in  the  prolonged  testing  needed  during  the.  combustor  design  and 
development  phase.  The  current  design  system  relies  on  considerable  empiricism  and 
"experience  database."  This  approach  has  limited  utility  when  radically  new  concepts  and 
designs  for  the  combustors  are  becoming  inevitable  due  to  the  increasingly  stringent 
requirements  and  aggressive  goals  for  the  combustors  as  well  as  the  engines  as  a  whole.  The 
limitations  of  the  current  design  system  stem  mainly  from  the  deficiencies  in  the 
physicochemical  models  for  the  key  combustion  processes  such  as  the  treatment  of  the  turbulent 
transport  of  momentum  and  scalars,  and  the  two-way  interactions  between  turbulence  and 
reactions.  Conventional  turbulent  combustion  models  suffer  from  the  turbulence  closure 
problem  (need  to  model  the  turbulent  transport)  and  the  inability  to  treat  realistic  finite-rate 
reactions  in  a  turbulent  flow  field.  A  fundamentally  different  approach  is  needed  to  address 
these  critical  issues. 

An  important  aspect  of  any  model  development  activity  is  the  validation  of  the  model  results 
against  reliable  benchmark  quality  data.  This  comparison  serves  to  validate  the  model  as  well  as 
to  guide  the  improvements  needed  for  the  model.  The  benchmark  experiments  must  be  designed 
and  carried  out  such  that  they  provide  accurate  inlet  and  boundary  conditions  as  well  as  the 
measurements  within  the  flow  field  for  all  the  quantities  needed  for  the  computations  and  model 
validation.  At  the  same  time,  the  combustor  geometries  should  be  simple  enough  to  permit  the 
accurate  simulation  of  the  flow  geometry  to  focus  attention  on  the  modeling  issues.  Finally,  the 
developed  model  needs  to  be  assessed  for  its  ability  to  compute  gas  turbine  combustor  flows. 
This  assessment  again  needs  to  performed  by  comparison  against  benchmark  quality 
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measurements  for  laboratory  combustors  which  incorporate  the  essential  features  of  gas  turbine 
combustors. 

The  current  program  seeks  to  address  the  above  issues  through  the  development  and  validation  of 
a  fundamentally  different  and  more  accurate  analytical  combustor  design  model.  The  program 
objectives  are  presented  in  the  next  section.  The  overview  of  the  program,  including  a  brief 
description  of  the  approach  and  the  key  achievements  under  the  program  are  given  in  the 
following  subsection.  An  overview  of  the  rest  of  this  report  is  also  presented  in  that  subsection. 


1.1  Program  Objectives 

The  objectives  of  this  program  are  the  following: 

•  Develop  and  validate  more  accurate  computer-based  analytical  models  for  combustor  flows. 

•  Design  benchmark  quality  experiments  and  use  the  data  from  advanced  diagnostic 
measurements  for  such  development  and  validation. 

•  Incorporate  the- analytical  models  into  a  next-generation  gas  turbine  combustor  design  system 
and  assess  its  capabilities. 


1.2  Overview  of  the  Present  Study 

The  joint  velocity-frequency-scalar  probability  density  function  (pdf)  method  (or  joint  pdf 
method  for  short)  is  the  focus  of  the  model  development  and  validation  efforts  at  Allison  under 
this  program.  The  pdf  method  was  selected  as  the  most  promising  candidate  for  further 
development  as  a  combustor  design  tool  after  a  review  of  several  current  and  evolving 
approaches.  The  pdf  method  overcomes  the  main  deficiencies  of  the  current  models;  turbulent 
transport  and  turbulence/chemistry  interactions  appear  in  closed  form  in  the  pdf  method  and  need 
not  be  modeled.  Before  the  start  of  the  current  program,  significant  progress  had  been  made  in 
the  development  of  the  pdf  approach  and  its  advantages  over  current  methods  had  been 
demonstrated  for  several  simple  reacting  and  nonreacting  flows.  However,  considerable  work 
needed  to  be  done  in  order  for  the  method  to  be  used  for  the  complex  swirling  recirculating  flow 
fields  in  practical  combustors. 

A  review  of  existing  benchmark  quality  data  showed  that  substantial  voids  existed  in  the  data 
with  respect  to  availability  of  inlet  condition  measurements  as  well  the  quality  and  completeness 
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of  the  data.  Also,  benchmark  quality  data  (meeting  the  criteria  described  above)  for  swirling  and 
recirculating  reacting  flows  were  almost  nonexistent.  Hence  benchmark  experimental 
configurations  were  designed  in  conjunction  with  Wright  Laboratory,  Combustion  Branch 
(WL/POSC)  personnel  and  University  of  Dayton  Research  Institute  (UDRI)  personnel. 
Measurements  were  made  by  UDRI  personnel  at  WL/POSC  under  a  separate  program.  These 
data  as  well  as  any  suitable  data  from  other  sources  (such  as  data  from  NASA  Hot  Section 
Technology  [HOST]  programs)  were  used  for  model  validation  in  the  current  program. 

Substantial  progress  has  been  made  under  the  current  program  in  developing  and  validating  the 
pdf  method  for  complex  flows  such  as  in  gas  turbine  combustors.  The  primary  areas  of 
progress  have  been  in  the  development  of  the  models  for  the  turbulent  frequency  and  the 
development  of  the  algorithms  for  determining  the  mean  pressure  field  within  the  pdf 
calculations.  Other  areas  have  included  the  treatment  of  swirling  flows  and  treatment  near  solid 
walls.  These  developments  are  key  to  the  application  of  the  pdf  method  to  complex  flows  of 
practical  interest.  The  pdf  method  has  been  applied  to  variety  of  swirling,  recirculating, 
nonreacting  and  reacting  flows  under  this  program.  Results  from  pdf  calculations  are  in 
excellent  agreement  with  detailed  benchmark  experimental  data  for  temperature  and  turbulent 
velocity  correlations  up  to  fourth  order.  The  results  from  the  joint  pdf  method  have  been 
assessed  against  the  results  from  other  conventional  methods  as  well  as  the  scalar  pdf  method  (a 
subset  of  the  joint  pdf  method).  The  pdf  methods  have  been  shown  to  be  more  accurate  than 
conventional  turbulent  combustion  models.  A  baseline  design  system  based  on  the  pdf  methods 
has  been  developed.  All  the  program  objectives  have  been  achieved. 

In  Chapter  2  of  this  report,  a  discussion  of  the  important  modeling  issues  for  gas  turbine 
combustor  flows  and  the  unique  advantages  offered  by  the  pdf  method  over  other  methods  is 
presented.  It  is  followed  by  a  description  of  the  pdf  method  including  the  governing  equations, 
the  modeling,  solution  algorithm  and  the  key  model  development  efforts  under  the  program.  The 
models  for  velocity  and  turbulent  frequency  developed  and  used  under  this  program  are 
described  in  Chapter  3.  The  models  used  for  the  molecular  mixing  of  scalars  are  described  in 
Chapter  4.  The  algorithms  developed  for  determining  the  mean  pressure  field  used  in  the  pdf 
solution  are  described  in  Chapter  5.  The  determination  of  the  mean  pressure  field  and  the 
associated  pdf  solution  algorithm  is  an  important  aspect  in  making  pdf  calculations  for  complex 
elliptic  recirculating  flows.  Chapters  6  through  10  describe  the  application  and  validation  of  the 
pdf  method  and  its  models  for  a  wide  range  of  flows  of  varying  complexities.  Available 
benchmark  quality  data  as  well  as  new  experiments  designed  and  conducted  to  provide 
benchmark  data  are  also  described  in  those  chapters.  The  joint  pdf  calculations  are  assessed 
against  results  form  the  scalar  pdf  method  and  several  conventional  models  in  Chapter  11. 
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Concluding  remarks  and  areas  of  further  development  of  the  pdf-based  design  system  are 
presented  in  Chapter  12. 
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2.  THE  PROBABILITY  DENSITY  FUNCTION  (PDF)  METHOD 


2.1  Background 

A  combustor  design  methodology  based  on  advanced  analytical  tools  can  be  used  with 
confidence  in  the  design  of  advanced  gas  turbine  combustion  systems  if  the  various 
physicochemical  processes  occurring  in  the  combustor  are  accurately  modeled  and  simulated. 
These  processes  include  the  mean  flow  dynamics;  turbulent  transport  of  momentum  and  scalars; 
molecular  mixing  of  scalars;  reaction  kinetics;  spray  injection,  droplet  formation,  dispersion,  and 
evaporation;  soot  formation  and  oxidation;  convective  and  radiative  heat  transfer;  and  the 
interactions  between  the  different  processes,  namely  turbulence/chemistry  interactions, 
turbulence/spray  dynamics  interactions,  etc.  In  addition,  accurate  numerical  schemes  and 
solution  algorithms  are  needed  to  solve  the  governing  equations  for  these  processes  to  obtain 
quantitatively  accurate  predictions  for  the  combustor  flow  and  performance  characteristics. 

The  model  development  and  validation  efforts  under  this  program  take  into  account  all  the  above 
aspects,  although  the  primary  attention  is  focussed  on  the  simulation  of  the  turbulent  reacting 
flow  field  (i.e.,  mean  flow,  turbulence,  scalar  transport,  reactions,  and  their  interactions).  Several 
past  and  ongoing  work  at  Allison  and  elsewhere  focus  on  the  advanced  spray  and  evaporation 
models  and  heat  transfer.  Attention  is  paid  in  the  development  of  the  method  in  the  current 
program  so  that  the  method  is  compatible  with  the  spray  and  other  advanced  models  and  that 
those  models  can  be  easily  incorporated. 

The  main  difficulties  in  computing  practical  turbulent  reacting  flows,  such  as  in  gas  turbines, 
arise  due  to  the  simultaneous  presence  of  turbulence,  mean  and  fluctuating  strain  rates  and 
pressure  gradients,  swirl,  body  forces,  scalar  transport,  complex  finite-rate  reactions  and  the 
interactions  between  the  various  processes.  Conventional  Reynolds-averaged  approaches,  which 
are  widely  used  currently,  suffer  from  the  "turbulence  closure  problem,"  inability  to  treat 
complex  finite-rate  reactions  and  to  accurately  account  for  the  interactions  between  turbulence 
and  reactions.  The  probability  density  function  (pdf)  method  used  and  developed  in  the  present 
study  are  remarkably  successful  in  alleviating  the  major  difficulties;  the  effects  of  turbulent 
transport,  complex  finite-rate  reactions,  body  forces,  variable-density,  mean  pressure  gradients, 
and  turbulence  chemistry  interactions  are  in  closed  form  and  can  be  treated  without 
approximation.  The  only  remaining  closure  problems  are  due  to  the  fluctuating  pressure  gradient 
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and  to  molecular  transport  (including  effects  of  fluctuating  velocity  and  scalar  gradients).  Some 
of  these  issues  and  the  advantages  of  the  pdf  method  are  discussed  in  more  detail  in  the  following 
subsections. 


2.1.1  The  Turbulence  Closure  Problem 

In  conventional  Reynolds-averaged  approaches,  turbulence  models  are  needed  to  overcome  the 
classical  turbulence  closure  problem,  i.e.  the  normal  and  shear  stresses  (the  so-called  Reynolds 
stresses)  are  unknown  in  the  mean  momentum  equations.  The  gradients  of  the  Reynolds  stresses 
that  appear  in  the  mean  momentum  equations  represent  the  turbulent  convective  transport  of 
momentum.  In  mean-flow  closure  models,  such  as  the  k-8  model  [2-1]*,  the  Reynolds  stresses 
are  related  to  the  mean  velocity  gradients  through  an  isotropic  eddy  viscosity,  based  on  the 
gradient  diffusion  (or  gradient  transport)  assumption.  Although,  robust  and  reasonably 
successful  in  predicting  a  variety  of  turbulence  flows,  the  k-8  model  has  serious  limitations  for 
accurately  predicting  the  location  and  extent  of  recirculation  zones  and  flows  with  significant 
streamline  curvature  such  as  in  curved  boundary-layers  and  swirling  flows,  since  the  assumption 
of  isotropic  eddy  viscosity  is  erroneous  for  such  flows  [2-2,  2-3].  Several  modifications  to  the  k- 
8  model  have  been  proposed  over  the  years  but  they  have  had  only  partial  successes  in 
overcoming  the  fundamental  problems  with  the  models. 

The  next  level  of  closure  in  the  Reynolds-averaged  approach  is  the  second-order  closure,  in 
which  transport  equations  for  the  Reynolds  stresses  are  solved  or  approximated  through  algebraic 
relations  [e.g.,  2-4,  2-5].  However,  in  these  equations  the  third-order  correlations,  which 
represent  the  turbulent  convective  transport  of  the  Reynolds  stresses,  are  unknown  and  need  to  be 
modeled.  Again,  the  modeling  approach  adopted  is  to  relate  the  third-order  correlations  to  the 
Reynolds  stresses  through  gradient  diffusion  assumptions.  The  second-order  closure  models 
better  account  for  the  anisotropy  of  turbulent  transport.  While  they  have  produced  better  results 
than  the  k-e  model  for  many  flows  especially  with  streamline  curvature,  they  have  not  been 
unequivocally  established  as  a  superior  alternative  to  the  k-e  model.  The  robustness  and  the 
simplicity  of  the  k-8  model  compared  to  the  second-order  models  make  the  former  a  more  widely 
used  choice  for  complex  practical  flows. 

The  main  drawback  of  both  classes  of  turbulence  models,  the  mean  and  second-order  closures,  is 
the  assumption  of  gradient  transport.  The  assumption  is  essentially  ad  hoc  and  is  not  valid  for 
many  flows,  especially  variable-density  and  reacting  flows  [2-6].  With  the  joint  velocity-scalar 
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pdf  method  [2-7],  the  turbulence  closure  problem  described  above  does  not  exist  and  hence  the 
gradient  transport  assumptions  are  avoided.  The  convective  transport  by  both  mean  and 
turbulent  velocities  appear  in  closed  form  and  need  not  be  modeled. 

An  alternative  approach  to  computing  turbulent  flows  is  direct  numerical  simulation  (DNS) 
[2-8].  DNS  attempts  to  simulate  the  stochastic  flow  field  exactly,  by  considering  the 
instantaneous  Navier- Stokes  equations.  To  truly  simulate  the  turbulent  processes,  the  numerical 
simulation  should  resolve  the  smallest  scales  of  turbulence.  This  fact  puts  severe  restrictions  on 
the  magnitude  of  the  Reynolds  number  of  the  flows  that  can  be  simulated  even  with  the  present 
day  high-speed  and  large-storage  computers.  DNS  cannot  be  used  for  computing  flows  in 
practical  geometries  of  interest,  much  less  as  a  design  tool  for  gas  turbine  combustors  for  many 
years  to  come.  The  utility  of  DNS  lies  largely  in  providing  guidance  for  developing  other 
turbulent  combustion  models.  Indeed  DNS  has  been  used  to  develop  the  turbulence  frequency 
models  developed  for  the  pdf  method  (see  Section  3.2).  In  another  approach,  large  eddy 
simulation  (LES),  only  the  large  energy-containing  scales  are  resolved  and  the  subgrid  scales  are 
modeled.  It  is  still  very  computer  intensive  and  impractical  for  design  calculations.  Moreover,  it 
has  the  same  limitations  as  conventional  approaches  for  calculating  reacting  flows  (see  next 
subsection). 


2.1.2  Chemistry/Reaction  Rates 

For  reacting  flows,  the  turbulence  models,  such  as  k-e  and  Reynolds  stress  models,  need  to  be 
extended  to  include  additional  models  for  scalar  flux,  mean  reaction  rates,  correlations  involving 
instantaneous  reaction  rates,  and  models  to  account  for  additional  effects  of  variable  density. 

Combustion  of  practical  fuels  typically  involve  several  hundred  coupled  reactions  and  tens  of 
species.  The  rates  of  these  reactions  vary  exponentially  with  temperature  and  are  non-linearly 
dependent  on  species  concentrations  (Arrhenius  rate  expressions).  In  Reynolds-averaged 
approaches  the  mean  reaction  rate,  i.e.,  the  mean  of  the  instantaneous  Arrhenius  reaction  rates 
based  on  the  instantaneous  temperature  and  composition  values,  is  needed  in  the  calculation. 
The  most  prevalent  current  practice  with  these  approaches  is  to  use  the  Arrhenius  reaction  rate 
corresponding  to  the  mean  temperature  and  composition  rather  than  the  mean  of  the 
instantaneous  rates.  Given  the  highly  nonlinear  reaction  rates  and  the  large  turbulent  fluctuations 
in  temperature  and  composition,  such  a  practice  can  lead  to  an  order-of-magnitude  error  in  the 
value  of  the  reaction  rate  and  can  even  yield  the  wrong  sign  [2-9].  The  current  combustion 
models  can  determine  mean  reaction  rates  accurately  only  under  very  restrictive  special  cases 
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(when  all  the  reaction  rates  are  linear  and  weakly  coupled  or  when  the  entire  combustion  process 
is  very  fast  or  very  slow).  A  general  nonlinear,  finite-rate,  multi-step  reaction  mechanism  cannot 
be  treated  by  these  conventional  models. 

The  most  striking  advantage  of  the  pdf  method  is  its  ability  to  treat  arbitrarily  complex, 
finite-rate,  coupled  multi-step  reactions  without  approximations.  (Although  this  is  theoretically 
true,  reasonable  simplifications  of  the  chemistry  may  be  required  to  reduce  computer  resource 
requirements.  There  is  ongoing  work  aimed  at  developing  efficient  computational  procedures 
without  requiring  simplification.) 

There  are  several  past  and  ongoing  efforts  to  overcome  the  reaction  closure  problem  in  the 
Reynolds-averaged  methods,  such  as  laminar  flamelet  model,  assumed  pdf  approaches, 
conditional  moment  closure,  etc.,  but  all  of  them  are  based  on  restrictive  simplifying  assumptions 
and  none  of  them  has  the  advantage  that  the  pdf  method  has  for  treating  reactions  in  the  turbulent 
combustion  of  practical  fuels  without  approximation.  Another  combustion  model,  the  linear 
eddy  model,  has  many  similarities  with  the  pdf  method  for  treating  the  reactions  but  lacks  the 
rigor  of  the  pdf  method  in  its  formulation  nor  is  it  directly  applicable  to  multidimensional  flows. 

2.1.3  Scalar  Transport 

In  addition  to  the  mean  reaction  rates,  the  calculation  of  scalar  transport  (i.e.,  transport  of  species 
and  other  scalars)  is  important  in  turbulent  reacting  flows.  In  conventional  Reynolds-averaged 
approaches,  the  turbulent  flux  of  the  scalar  needs  to  be  modeled  while  solving  the  conservation 
equation  for  the  scalar.  The  same  ideas  behind  the  modeling  of  the  Reynolds  stresses  (section 
2.1.1),  namely  gradient  diffusion,  are  applied  for  the  modeling  of  scalar  fluxes  as  well.  As  with 
the  momentum  flux,  the  gradient  diffusion  assumptions  for  scalar  flux  can  be  erroneous 
especially  in  reacting  flows  (see  section  2.1.4).  With  the  joint  pdf  method,  the  (mean  and) 
turbulent  flux  of  scalar  is  in  closed  form  and  need  not  be  modeled. 

2.1.4  Turbulence/Chemistry  Interactions 

Several  processes  such  as  turbulence  production  and  dissipation,  turbulent  transport,  reaction  and 
the  accompanying  heat  release  and  volume  expansion,  occurring  in  turbulent  combustion  are 
strongly  interconnected.  The  effect  of  turbulent  fluctuations  on  mean  reaction  rates  was 
described  in  Section  2.1.2.  In  turn,  the  turbulence  is  influenced  by  the  reaction. 

The  volume  expansion  (density  change)  resulting  from  reaction  affects  turbulence  production 
and  dissipation  through  changes  in  velocity  gradients  and  interactions  between  variable-density 
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and  the  mean  pressure  gradients.  While  these  effects  are  difficult  to  model  in  the  conventional 
approaches,  they  appear  in  closed  form  in  the  joint  pdf  method  and  need  not  be  modeled. 

Another  important  effect  of  the  heat  release  is  the  occurrence  of  the  phenomenon  of 
counter-gradient  diffusion  [2-10,  2-11]  in  which  the  scalar  flux  and  the  scalar  gradient  are  of  the 
same  sign  in  certain  regions  of  the  flow  contrary  to  conventional  gradient  transport  modeling 
assumption.  A  more  general  case  of  this  phenomenon  is  when  the  scalar  flux  vector  is  aligned  at 
an  angle  other  than  diametrically  opposite  to  the  scalar  gradient  [2-12,  2-13].  The  reason  for  this 
phenomenon  lies  in  the  fact  that  the  difference  between  the  densities  of  the  reactants  and  the 
products  can  be  quite  large,  and  the  strong  mean  pressure  gradient  across  the  reaction  zone 
imparts  greater  acceleration  to  the  low-density  products  than  to  the  reactants.  None  of  the 
conventional  models  can  treat  this  phenomenon  in  a  general  way.  Special  models  are  used  for 
particular  situations  under  consideration  [2-10,  2-14].  With  the  joint  pdf  method,  the 
phenomenon  of  counter-gradient  diffusion  and  other  variable-density  effects  are  automatically 
calculated  as  demonstrated  by  Anand  et  al.  [2-11, 2-13]. 

In  summary,  the  joint  pdf  method  is  inherently  successful  in  overcoming  the  most  important 
modeling  difficulties,  such  as  turbulent  transport,  complex  reactions,  and  turbulent  chemistry 
interaction.  The  pdf  method  provides  a  more  complete  description  of  the  turbulent  flow  field 
than  the  conventional  models  that  usually  calculate  only  the  first  and  second  moments.  In 
addition,  when  there  are  multiple  streams  or  fluids  conditional  modeling  is  easily  accomplished 
in  the  pdf  method.  Effects  of  large  scale  structures  and  phenomena  such  as  intermittency, 
unmixedness,  etc.  can  be  accounted  for.  These  advantages  make  the  pdf  method  the  most 
suitable  approach  for  computing  turbulent  reacting  flows  and  warrant  its  development  as  a  gas 
turbine  combustor  design  tool. 


2.2  Description  of  the  PDF  Method 

There  is  a  hierarchy  of  pdf  methods  currently  in  use.  They  are  based  on  (a)  the  joint  pdf  of 
compositions  (b)  the  joint  pdf  of  velocity  and  composition,  U,  &  and  (c)  the  joint  pdf  of 
velocity,  turbulent  frequency  and  composition,  U,  co,  <|).  Methods  (b)  and  (c)  are  referred  to  as 
the  joint  pdf  method  in  this  report  and  method  a)  is  referred  to  as  the  scalar  pdf  method.  The 
primary  focus  of  the  present  study  is  on  the  most  advanced  and  comprehensive  method  (c),  while 
work  with  method  (b)  has  also  been  performed.  The  scalar  PDF  method  overcomes  the  difficulty 
with  treating  complex  chemistry/reaction  rates  (section  2.1.2)  and  the  effect  of  turbulence  on 
chemistry.  The  joint  pdf  method  in  addition  overcomes  the  turbulent  transport  issue  (of 
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momentum  and  scalars)  and  more  comprehensively  treats  the  two-way  interactions  between 
turbulence  and  chemistry  (section  2.1.4). 

The  joint  pdf  approach  involves  the  solution  of  the  transport  equation  of  the  velocity-frequency- 
scalar  pdf.  The  joint  PDF  f(V,  T|,  vjr;  x,  t),  is  the  probability  density  of  the  occurrence  of  the 
simultaneous  event  U(  x,  t)  =  V,  (0  (x,  t)  =  T|  and  £(x,  t)  =  \g,  where  U  is  the  velocity  vector,  co  is 
the  turbulence  frequency,  $  is  a  set  of  scalars  and  V,  r\,  and  \j£  are  independent  variables  in  the 
velocity,  frequency  and  scalar  spaces  respectively.  Any  single  point  correlation  (means, 
variances,  third  moments,  etc.)  at  a  position  x  and  time  t  can  be  evaluated  from  the  following 
properties  of  the  joint  pdf: 


Q(U,0)^)  =  JJJ  Q(V,  T),  \jr)  f(V,  r),  \jt;  x,  t)  dVdTid\ir, 

(2.1) 

JJJ  f(Y,  ti,  \jr;  x,  t)  dV  dr]  d\j£  =  1  • 

(2.2) 

The  integrals  in  Equations  2.1  and  2.2  are  over  the  entire  velocity-frequency- scalar  space;  Q 
represents  any  quantity  (or  correlation)  whose  mean  or  "expectation"  is  required;  and  the  overbar 
in  Equation  2. 1  represents  the  Reynolds-averaged  mean. 


For  variable-density  reacting  flows,  it  is  more  appropriate  to  consider  the  density-weighted  pdf  f 
or  alternatively  the  mass  density  function  (mdf)  F  defined  by 

f(V,  ti,  \j/;  x,  t)  =  pOj/)  f(V,  ti,  n/;  x,  t)  /  p  ,  (2.3) 


and 

F(V,  T|,  \jr;  x,  t)  =  p(\j£)  f(V,  T],  \jr;  x*  t)  =  p  f  ,  (2.4) 

respectively,  where  p  is  the  mean  density.  The  density-weighted  or  Favre-averaged  mean 
(denoted  by  angled  brackets)  of  any  quantity  Q  is  given  by: 

<Q(U,co,<!>)>  =PQ/P 

=  JJJ  P(¥)QCV,ti,¥)  f(Y,ti,3|£)  dV  dr]  d\|/  /  p 
=  JJJ  Q(Y,  ti,  \jr)  f(V,  ti,  \j/)  dV  dr\  d\j£ 
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=  JJJ  Q(Y,  T|,  \jr)  F(V,  T|,  \jr)  d V  dr)  d\jr  /  p  . 


(2.5) 


(Note:  The  attributes  x  and  t  are  implied  and  not  specifically  stated  for  brevity.) 

The  pdf  f  represents  the  number  probability  density  of  events  in  the  state  space  (velocity- 
frequency-scalar  space)  while  the  mdf  F  represents  the  probability  density  of  the  mass  associated 
with  the  events. 


2.2.1  The  PDF  Transport  Equation 

The  transport  equations  for  f,  f  or  F  can  be  derived  starting  from  the  familiar 
instantaneous  conservation  equations  for  mass,  momentum,  scalar  (species,  enthalpy,  etc.)  and  a 
general  model  equation  for  the  turbulence  frequency.  The  conservation  equations  are  listed 
below. 


Continuity  (instantaneous): 


9(pUj) 

3t +  dxj 


=  0  , 


(2.6) 


where  p(^)  is  the  instantaneous  fluid  density  uniquely  determined  by  the  set  of  scalars  <j),  and 
U[  is  the  instantaneous  velocity  in  the  i-th  direction. 

Momentum  (instantaneous): 


DUj  1  9P 


Dt  p  9xj 


^7  +gi  +Ai  ’ 


where  the  material  derivative  (or  the  rate  of  change)  following  the  fluid  is 


D _ d_  d_ 

Dt  3t  +  k  dx. 


(2.7) 


(2.8) 


The  first  term  on  the  right-hand  side  (of  Eq.  2.7)  represents  the  mean  pressure  gradient,  gj  is 
the  body  force  per  unit  mass  in  the  j-th  direction  and 
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J  p  W*.  dx}J 


(2.9) 


represents  the  terms  due  to  the  sum  of  the  viscous  and  viscous-diffusive  stress  tensor  (xy) 
and  the  fluctuating  pressure  gradient  (BpVdxj).  Note  that  the  instantaneous  pressure  P  has 
been  decomposed  into  its  mean  and  fluctuating  parts  while  writing  Eq.  2.7,  i.e. 


P  =  P+p'  . 


(2.10) 


Scalar  conservation  (instantaneous): 


Dt 


Sa  +  0a  *  cc  —  1,  2, 


c y , 


(2.11) 


where  the  number  of  scalars  considered  is  o,  Sa  is  the  rate  of  addition  of  scalar  a  (per  unit 
mass  of  the  fluid)  due  to  reaction  or  other  sources,  and 


(2.12) 


where  Ja  is  the  diffusive  flux  vector  of  species  a.  Note  that  the  scalars  could  represent 
species  mass  fractions,  enthalpy,  conserved  scalars  such  as  mixture  fraction,  etc. 

Turbulence  frequency  (instantaneous): 


Deo 

Dt 


A  , 


(2.13) 


where  A,  which  can  be  derived  from  the  Navier- Stokes  equation,  describes  the  evolution  of 
the  turbulent  frequency.  The  turbulence  frequency  is  the  inverse  of  the  turbulence  time  scale 
x  and  is  defined  by 


co  =  x-1  =  e/k 


(2.14) 


where  k  is  the  turbulent  kinetic  energy  and  e  is  the  rate  of  dissipation  of  k. 
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The  molecular  transport  quantities  and  _Ja,  and  A  are  complicated  functions  of  the  local 
properties  and  their  gradients. 

The  transport  equation  for  the  mdf,  derived  from  the  conservation  equations  2.6,  2.7,  2.11,  and 
2.13  using  simple  techniques  described  in  Ref.  2-7,  is: 


0F  l  1  d  [o  , 

at  ^  +(gj  Pty)  axj/avj  J 


av; 


a¥a' 


9  [<Aj|y,Ti,tjr>F]  ~ [<©a|V,T1,t!/>F]  -|j-[<A|V,tl,xir>F]  . 


(2.15) 


The  terms  on  the  left-hand  side  of  the  mdf  transport  equation  represent,  respectively,  the  time 
rate  of  change  of  the  mdf,  transport  in  physical  space  (convection),  transport  in  velocity  space  by 
mean  pressure  gradient  and  body  forces,  and  transport  in  scalar  space  due  to  reaction  for 
example.  All  the  terms  on  the  left-hand  side  are  known  functions  of  the  independent  variables  x, 
t,  V,  r),  and  \jt.  Hence  they  are  in  closed  form  and  need  not  be  modeled.  The  terms  on  the 
right-hand  side  involve  conditional  expectations.  For  example,  <Aj  I V,  T|,  \j£>  is  the  conditional 
expectation  of  Aj  given  that  U  =  V,  to  =  i),  and  ^  The  conditional  expectations  involve  the 
effects  of  viscous  dissipation,  molecular  diffusion  and  the  fluctuating  pressure  gradient.  They 
are  in  general  unknown  and  need  to  be  modeled.  It  should  be  noted  that  these  effects  are 
modeled  in  the  conventional  methods  as  well.  Hence  there  is  no  additional  modeling  burden 
created.  In  fact,  the  major  modeling  difficulties  with  the  conventional  methods  are  alleviated  in 
the  pdf  method.  The  pdf  method  also  allows  a  more  sophisticated  and  detailed  modeling  of  the 
molecular  and  fluctuating  pressure  gradient  terms  than  conventional  models. 


2.2.2  Modeling  and  Solution  of  the  PDF  Transport  Equation 

.  The  transport  equation  for  the  pdf  (or  the  mdf)  could  in  principle  be  solved  using  finite- volume 
methods.  However,  due  to  the  large  dimensionality  of  the  pdf  (a  total  of  7+cj  independent 
variables  including  spatial  coordinates,  time,  velocity  and  scalars)  finite-volume  methods  are 
uneconomical  since  the  computer  requirements  rise  exponentially  with  the  dimensionality  of  the 
pdf.  Monte  Carlo  methods  are  most  suited  and  economical  for  these  computations.  In  the  Monte 
Carlo  method  used,  the  flow  field  is  simulated  using  a  large  number  of  computational  particles 
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(representing  fluid  particles)  whose  evolution  mimics  the  evolution  of  fluid  particles.  Hence  a 
Lagrangian  viewpoint  is  adopted  in  modeling  and  solving  the  mdf  transport  equation. 

In  the  solution  technique,  notional  particles  representing  fluid  particles  are  distributed  throughout 
the  solution  domain  overlaid  by  a  spatial  or  computational  grid.  Each  particle  represents  a 
specified  fluid  mass  and  is  attributed  with  values  for  its  spatial  position  (x*),  velocity  (U*), 
scalar  values  ($*)  and  turbulence  frequency  (co*).  These  values  evolve  according  to  the 
equations  described  below  which  include  modeled  terms  where  needed.  Starting  from  arbitrary 
initial  conditions  and  specified  boundary  conditions,  the  particle  values  are  marched  in  time- 
steps  which  are  a  fraction  of  a  characteristic  time  scale  in  the  flow  until  a  steady-state  solution  is 
reached.  As  will  be  demonstrated  below,  the  solution  of  these  evolution  equations  constitutes 
the  solution  of  the  mdf  transport  equation.  Means  of  any  functions  of  the  independent  variables 
are  determined  by  summing  the  values  of  the  required  properties  over  all  the  particles  in  each  of 
the  spatial  bins  or  cells  and  fitting  splines  through  these  sums  over  the  whole  domain.  More 
details  can  be  obtained  from  Ref.  2-7  and  other  works  listed  in  the  Appendix. 

The  increment  dx*  in  the  position  of  a  particle  over  an  infinitesimal  time  interval  dt  during  the 
time  step  is  given  by  the  exact  equation: 


dx*  =  U*  dt  .  (2-16) 

This  exact  equation  causes  the  mean  and  turbulent  convection  to  be  in  closed  form  since  the 
velocity  used  is  the  instantaneous  velocity  which  includes  the  mean  and  fluctuation. 

The  general  form  for  the  increment  in  the  particle  velocity  is 

dU*  =  -  dt  +gj  dt  +a;  dt  ,  (2.17) 

i  p*  dx.  1 

where  p*($*)  is  the  density  of  the  particle,  and  a  represents  the  model  for  the  effects  of  the 
viscous  dissipation  and  fluctuating  pressure  gradient.  In  general,  the  model  a  consists  of  both 
deterministic  and  random  parts.  The  model  a  is  constructed  such  that  it  is  determined  by  the 
current  properties  of  the  computational  particles.  In  other  words,  a  is  fully  evaluated  from  the 
current  mdf  but  includes  a  random  component. 

Similarly,  the  general  forms  of  the  increments  in  the  scalar  and  frequency  values  of  the  particle 
over  an  infinitesimal  time  interval  dt  are: 
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#«*  =  Sa(<|>*)  dt  +  0a  dt , 


(2.18) 


and 

dco*  =  Adt  .  (2.19) 

Again,  the  models  fi  and  A  could  have  deterministic  and  random  parts  and  are  constructed  such 
that  they  are  evaluated  from  the  current  mdf. 

The  transport  equation  of  the  mdf  F  (which  is  an  Eulerian  quantity  by  definition)  can  be  derived 
from  the  Lagrangian  particle  evolution  equations  (2.16  through  2.19)  using  techniques  presented 
in  Ref  2-7.  The  resulting  equation  is: 

8F  3F  /  1  BP\  BF  B  r_  , 

at  +  VjdXj  +  lgj  P(\jt)  BxJ  BVj  +  -* 

-  5^[ajlI’Tl’21f>F]  -a“[<0«IXTl,}l/>F]  -^[<^|V,Ti,i|/>F]  .  (2.20) 

A  comparison  of  Eqs.  2.15  and  2.20  shows  that  the  mdf  of  the  computational  particles  (Eq.  2.20) 
evolves  identically  as  the  mdf  of  the  fluid  provided  the  corresponding  conditional  expectations 
(i.e.,  those  of  a  and  A,  fi  and  £),  and  A  and  A)  are  equal.  In  other  words,  as  long  as  the  models 
are  constructed  such  that  "on  the  average"  the  computational  particles  behave  the  same  way  as 
the  fluid  particles,  the  solution  of  the  evolution  equations  (2-16  through  2-19)  is  statistically 
equivalent  to  solving  the  mdf  transport  equation  for  the  fluid  and  consequently  the  conservation 
equations  for  the  fluid.  Any  turbulent  correlation  of  any  order  extracted  from  the  computed  mdf 
will  accurately  reflect  the  correlations  in  the  fluid  even  though  the  computational  particles  may 
behave  differently  from  the  fluid  particle  in  microscopic  detail. 

Some  of  the  models  for  a  and  A  are  described  in  the  next  chapter  and  the  models  for  Q  are 
described  in  chapter  4.  These  models  typically  need  information  of  the  turbulence  time  scales. 
The  motivation  for  including  the  turbulence  frequency  and  the  development  of  the  joint 
velocity-frequency-scalar  pdf  method,  is  that  the  turbulence  time  scales  are  determined  within 
the  calculations. 
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2.3  Previous  Work 


A  review  of  the  development  and  application  of  the  pdf  method  prior  to  the  present  program  can 
be  found  in  Refs.  2-7,  2-13  and  2-15.  Although  there  had  been  some  preliminary  work  with  the 
velocity  pdf  in  the  late  sixties  by  Lundgren  and  some  scalar  pdf  work  in  the  mid-seventies  by 
Dopazo  and  O’Brien,  and  Pope,  the  development  of  the  joint  velocity-scalar  pdf  method  was 
pioneered  by  Pope  in  the  early  eighties.  Pope  and  coworkers,  including  Anand  et.  al.,  applied  the 
joint  pdf  method  to  several  simple  reacting  and  nonreacting  flows.  The  flows  studied  included 
homogeneous  turbulence,  self-similar  jets,  wakes,  and  mixing  layer,  the  thermal  wake, 
one-dimensional  planar  and  spherical  premixed  flames,  and  diffusion  flames.  All  these  flows 
were  either  two-dimensional  (2-D)  parabolic  or  1-D  time-dependent  flows.  These  fundamental 
studies  contributed  to  the  systematic  and  rigorous  development  of  the  pdf  method,  including  its 
basic  models  and  solution  algorithms,  and  demonstrated  the  advantages  (discussed  in  section  2.1) 
of  the  pdf  method  for  computing  turbulent  reacting  flows. 

All  previous  work  involved  at  most  the  sophistication  of  the  joint  velocity-scalar  pdf  method.  A 
major  drawback  of  the  joint  velocity-scalar  pdf  method,  i.e.,  without  the  turbulence  frequency,  is 
that  it  does  not  contain  any  information  on  the  time  scales  of  the  turbulent  motions,  needed  for 
the  models  used  in  the  calculations.  However,  given  the  time  scale  information,  the  method 
provides  a  complete  description  of  the  one-point  joint  velocity-scalar  statistics  throughout  the 
flowfield.  Previous  studies  discussed  in  the  previous  paragraph  assumed  either  a  uniform 
turbulence  time  scale  across  the  flow  that  was  estimated  from  mean  flow  properties  or  otherwise 
externally  supplied  the  turbulence  time  scale  to  the  pdf  calculations.  The  assumption  of  a 
uniform  time  scale  is  not  valid  in  general  except  in  homogeneous  turbulence  or  in  specifically 
designed  experiments.  In  fact,  in  most  practical  reacting  and  non-reacting  flows  the  turbulence 
time  scale  varies  significantly  within  the  flow. 

Hence  a  major  thrust  of  the  present  study  has  been  toward  the  development  and  validation  of  the 
joint  velocity-frequency-scalar  pdf  method  in  which  the  turbulence  frequency  (inverse  of  the 
turbulence  time  scale,  which  is  the  primary  quantity  of  interest)  is  also  included  as  a  random 
variable  in  the  joint  pdf.  This  causes  the  joint  pdf  equation  and  the  calculations  to  be  fully 
self-contained. 

Another  important  issue  in  the  joint  pdf  method  is  the  calculation  of  the  mean  pressure  gradient 
needed  for  solving  the  evolution  equation  for  the  particle  velocities.  Theoretically,  there  is  no 
problem  in  determining  the  mean  pressure  gradients.  Within  every  iteration  or  time-step  in  the 
calculations,  the  mean  pressure  field  can  be  extracted  from  the  pdf  through  the  solution  of  the 
Poisson  equation  for  pressure.  However,  a  robust  algorithm  is  needed  to  solve  the  Poisson 
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equation  which  involves  the  evaluation  of  second  derivatives  of  mean  velocities  and  other  terms 
which  are  very  sensitive  to  statistical  errors  inherent  in  the  calculations. 

For  parabolic  flows,  the  boundary-layer  assumptions  can  be  used  to  determine  the  transverse 
pressure  gradient  from  the  transverse  momentum  equation  (see  section  5.1).  A  pdf  method  for 
elliptic  flows  was  first  demonstrated  by  Anand  et  al.  [2-15]  for  a  2-D  flow.  The  method  was 
subsequently  applied  to  3-D  in-cylinder  flows  for  spark  ignition  engines  by  Haworth  and  El 
Tahry  [2-16, 2-17].  In  this  method,  the  Monte  Carlo  (MC)  calculations  for  the  pdf  are  combined 
with  conventional  finite-volume  (FV)  calculations  [2-18]  for  the  mean  fields.  The  mean  pressure 
field  and  the  turbulence  time  scale  are  supplied  to  the  MC  calculations  by  the  FV  calculations. 
In  turn,  the  MC  calculations  supply  to  the  FV  calculations  the  Reynolds  stresses  extracted  from 
the  pdf  solution,  thereby  avoiding  the  use  of  conventional  turbulence  models.  This  combined 
algorithm  exploits  the  advantages  of  both  methods  to  yield  an  economical  algorithm  for  pdf 
calculations  of  elliptic  flows.  This  method  was  intended  to  demonstrate  the  feasibility  of  pdf 
calculations  for  such  flows. 

However,  for  more  complex  flows,  especially  variable  density  and  reacting  flows  with  fast  or 
finite-rate  multistep  chemistry,  the  coupling  between  the  two  methods  becomes  quite  complex 
and  it  becomes  more  difficult  to  fully  exploit  the  advantages  of  the  pdf  method.  Therefore,  it  is 
desirable  to  perform  pdf  calculations  independently,  without  recourse  to  FV  calculations.  Hence, 
another  major  aspect  of  the  present  study  is  to  develop  and  validate  a  robust  and  efficient 
algorithm  for  determining  the  mean  pressure. 

A  third  aspect  of  the  present  study  that  builds  on  the  previous  studies  is  the  application  of  the  pdf 
method  to  more  complex  flows  of  practical  interest  and  validation  against  experimental  data.  As 
explained  previously,  the  previous  studies  were  involved  with  fundamental  development  of  the 
pdf  method  and  models  and  served  to  demonstrate  the  advantages  of  the  pdf  method. 
Understandably,  the  studies  involved  many  simple  and  idealized  flows.  Results  were  compared 
with  data  from  fundamental  experiments  and/or  known  behavior  in  flows  such  as  homogeneous 
turbulence,  decaying  grid-generated  turbulence  self-similar  flows  and  some  simple  flames.  The 
method  needed  to  be  applied  to  and  validated  for  more  complex  flows  involving  the  essential 
features  of  gas  turbine  combustor  flows  such  as  swirl,  recirculation  and  realistic  multi-step 
chemistry,  for  the  method  to  be  developed  as  a  gas  turbine  computer  design  tool. 
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2.4  Key  Accomplishments  under  the  Present  Study 


The  objective  of  the  present  program  was  to  perform  the  necessary  model  and  algorithm 
development  and  validation  to  develop  the  pdf  method  as  a  combustor  design  and  analysis  tool. 
Consistent  with  the  discussion  in  the  previous  subsection,  two  of  the  major  aspects  of  the  present 
study  are  the  development  of  the  joint  velocity-frequency-scalar  method,  and  the  mean  pressure 
algorithms  for  elliptic  flows.  Three  variants  of  the  frequency  and  velocity  models  and  two 
variants  of  the  pressure  algorithm  were  developed  and  validated.  A  method  for  computing 
swirling  flows  with  the  boundary-layer  algorithm  was  developed  and  validated.  The  validation 
of  the  models  were  performed  for  a  wide  variety  of  flows  starting  with  fundamental  flows  for 
model  development  and  preliminary  validation  as  well  as  more  complex  gas-turbine-combustor¬ 
like  flows  for  final  validation.  Available  benchmark  quality  data  were  used  where  possible. 
However,  due  to  incompleteness  of  data  sets  or  lack  of  good  quality  data  for  flows  of  practical 
interest,  benchmark  quality  experiments  were  designed  in  conjunction  with  WPAFB  and  the 
University  of  Dayton  Research  Institute  (UDRI)  personnel.  The  experiments  were  set  up  and 
performed  by  UDRI  personnel  at  WPAFB. 

The  specific  key  accomplishments  are: 

•  development  of  the  model  for  mean  turbulent  frequency  and  validation  for  general 
inhomogeneous  turbulent  flows 

•  development  of  the  treatment  of  swirl  in  the  boundary-layer  algorithm 

•  development  of  the  stochastic  models  for  turbulence  frequency  and  the  corresponding  models 
for  the  velocity 

•  validation  of  the  above  for  a  variety  of  flows,  including  non-swirling  and  swirling, 
non-reacting  and  reacting  flows,  through  comparison  against  benchmark  quality  data 

•  development  and  validation  of  the  treatment  near  walls  for  wall  bounded  flows 

•  development  and  validation  of  the  mean  pressure  algorithm  for  elliptic  flows 

•  development  of  an  improved  mean  pressure  algorithm 

•  design  of  benchmark  experiments.  (The  data  are  not  specific  to  the  pdf  method.  They  can  be 
used  by  the  scientific  community  for  validation  of  other  models  as  well.) 

•  validation  against  backward-facing  step  flow  and  step-swirl  methane  combustor  (swirling, 
recirculating  and  reacting  flows) 

•  comparison  of  joint  pdf  calculations  against  calculations  from  scalar  pdf  and  other 
conventional  models  (design  system  assessment  performed  under  Allison  IR&D  program) 

•  publication  of  several  papers.  (A  list  of  the  papers  is  presented  in  the  Appendix.) 
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The  pdf  model  development  and  validation  efforts  have  been  very  successful.  The  results  of  pdf 
calculations  are  in  very  good  agreement  with  data  for  temperature  and  velocities  (mean,  variance 
and  up  to  fourth  order  turbulent  correlations).  A  baseline  version  of  the  joint  pdf  based 
combustor  design  system  has  been  developed.  Several  enhancements  to  the  design  system  have 
been  identified  and  planned  for  the  future  to  fully  meet  the  requirements  of  the  gas  turbine 
combustor  design  system. 
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3.  MODELS  FOR  VELOCITY  AND  TURBULENCE  FREQUENCY 


this  chapter  describes  the  models  developed  and/or  used  in  the  present  program,  for  velocity  and 
turbulence  frequency  in  the  particle  evolution  equations  (2.17  and  2.19  respectively).  Three 
different  models  for  turbulence  frequency  and  the  corresponding  velocity  models  are  described 
and  discussed  in  the  sections  to  follow. 


3.1  Mean  Frequency  Model  (MFM) 

The  mean  frequency  model  was  developed  as  a  means  for  supplying  the  turbulence  time  scale  to 
the  joint  velocity-scalar  pdf  calculations  before  the  joint  velocity-frequency-scalar  pdf  method 
was  developed.  In  this  model  the  randomness  of  the  turbulence  frequency  go  is  not  accounted 
for;  rather  a  transport  equation  for  the  mean  turbulence  frequency  <co>  is  solved  in  conjunction 
with  the  pdf  calculations  to  supply  the  time  scale  for  the  pdf  calculations.  Details  of  the 
development  of  the  model  are  presented  by  Anand  et  al.  [3-1]. 


The  transport  equation  for  <co>  (=  k/<e>)  is  derived  from  the  standard  equations  [3-2]  for  the 
turbulent  kinetic  energy  k  and  the  mean  dissipation  <e>  with  a  simplification  of  the  resulting 
diffusion  terms: 


0<to> 

at 


+  <Ui>^  =  CjSjjSjj  -  C2<co>2  +  -* 
9x;  1  ,J  1J  z 


a 

k  a<co>- 

ax- 

.<©>  ax;  . 

(3.1) 


where  Sy  is  the  mean  rate  of  the  strain  given  by 

1  /  a<U:>  a<U:>  \ 

Sij  =  2V  dXj  +  ax,  )  ’ 


(3.2) 


and  Ci  (=0.1),  C2  (=0.88),  (=0.09)  and  ow  (=1.0)  are  constants.  The  values  of  the  constants 

are  derived  from  the  standard  values  of  the  constants  in  the  k  and  <e>  equations.  The  left-hand 
side  of  Eq.  3.1  represents  the  transport  of  <co>  by  the  mean  flow.  The  terms  on  the  right 
represent  the  production,  dissipation  and  turbulent  diffusion  of  <0)>  respectively.  Eq.  3.1  is 
solved  by  finite  differences  using  the  mean  velocity  and  kinetic  energy  fields  extracted  from  the 
pdf  solution.  The  mean  frequency  in  turn  is  used  in  the  models  used  in  the  pdf  method.  This 
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frequency  model  was  used  in  conjunction  with  the  boundary  layer  algorithm  (section  5.1)  and 
applied  for  the  calculations  presented  in  Chapters  6  and  7. 


3.2  Simple  Langevin  Model  I  for  Velocity  using  MFM  (SLM I) 

The  velocity  model  used  belongs  to  a  class  of  stochastic  models  called  the  Langevin  model.  The 
models  used  here  are  a  simplification  of  the  generalized  Langevin  models  used  in  some  previous 
studies  with  the  pdf  method  [3-3].  The  simple  Langevin  model  used  by  Anand  and  Pope  [3-4]  is 
considerably  simpler  to  implement  and  sufficiently  accurate  for  most  of  the  flows  under 
consideration. 

The  simple  Langevin  model  for  velocity  corresponding  to  the  mean  frequency  model  (section 
3.1)  is: 

dU*  =  -  dt  +g,  dt  -  (|+|c<J)(u*-<Ui>)<o»dt  +(C,1<d»k)1'2dWi  ,  (3.3) 

where  C0  is  a  universal  constant  (due  to  Kolmogorov),  and  dWj  represents  an  isotropic  Wiener 
random  process  with  the  properties  that  its  mean  and  variance  have  the  values  <dWj>  =  0  and 
<dWjdWj>  =  dt  5jj.  A  comparison  of  Eq.  3.3  with  Eq.  2.17  shows  that  the  last  two  terms  in  Eq. 

3.3  correspond  to  a;  in  Eq.  in  2.17  and  collectively  model  the  effects  of  viscous  dissipation  and 
fluctuating  pressure  gradients.  The  first  of  the  two  terms  is  a  deterministic  term,  called  the  drift 
term  since  it  tends  to  relax  the  particle  velocity  toward  the  mean  velocity  and  reduce  the  velocity 
variance.  The  second  of  the  two  terms  is  a  random  term  and  models  the  retum-to-isotropy 
phenomenon  of  turbulence.  The  constant  CD  was  determined  by  Anand  and  Pope  [3-4]  to  have 
the  value  2.1.  As  discussed  in  section  2.2.2  the  first  two  terms  on  the  right-hand  side  of  Eq.  3.3 
are  exact  and  account  for  the  effects  of  the  mean  pressure  gradient  and  body  forces  without 
requiring  modeling.  Note  that  the  first  term  contains  the  interaction  between  the  instantaneous 
particle  density  and  the  mean  pressure  gradient.  This  is  the  reason  that  turbulence  chemistry 
interaction  phenomena  such  as  counter-gradient  diffusion  and  turbulence  modulation  due  to 
reaction  (see  section  2.1.4)  are  automatically  calculated  in  the  joint  pdf  method. 


3.3  Stochastic  Frequency  Model  I  (SFM I) 

The  stochastic  frequency  model  provides  the  foundation  for  the  joint  velocity-frequency-scalar 
pdf  method.  The  turbulence  frequency  is  considered  a  random  property  in  the  calculations  which 
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is  consistent  with  the  real  situation.  Indeed  turbulent  flow  consists  of  a  range  of  length  and  time 
scales.  In  addition,  the  model  accounts  for  internal  intermittency  of  turbulence,  accounts  for  the 
effects  of  large  scale  turbulent  structures,  and  better  accounts  for  the  origin  and  history  of  the 
fluid  particles. 

The  stochastic  frequency  model  was  constructed  with  reference  to  known  and  computed  (by 
direct  numerical  simulations  [3-5])  properties  of  homogeneous  turbulence.  The  model  was  then 
extended  to  inhomogeneous  turbulence  and  validated  against  several  self-similar  flows.  The 
model  was  then  extensively  exercised  and  validated  for  more  complex  practical  flows.  A  brief 
description  of  the  model  is  presented  here.  Details  of  the  construction  and  validation  of  the 
model  are  presented  in  Refs.  3-6  through  3-9. 

According  to  the  model,  the  increment  in  the  particle  frequency  over  an  infinitesimal  time 
interval  dt  is  given  by 

dm*  =  7“CcoiSijSijdt  “  ®*<©>[Cffl2+CxQ/n]dt  +co*(2Cz<o)>o2)1/2dW  +<©>2hdt, 

(3.4) 

where 

Q[n  =  ln(  (£>*/<(£»)  -  <  (co*/<co>)  ln( to*/<co>)  >  ,  (3.5) 

and  dW  represents  a  Wiener  random  process.  The  constants  in  the  equation  (3.4)  have  the  values 
=  1.6  and  a  =  1.0.  The  values  Cwi  =  0.04  and  Ca2  =  0-70  were  found  to  result  in  the  best 
agreement  with  data  (see  Chapter  7).  The  term  involving  the  quantity  h  serves  to  entrain 
nonturbulent  fluid  in  the  intermittent  region  and  is  discussed  in  detail  in  Ref.  3-7.  Additional 
comments  regarding  the  role  of  the  term  are  presented  in  Chapter  7.  This  model  has  been  used 
for  computing  flows  presented  in  Chapters  7  and  8. 


3.4  Refined  Langevin  Model  for  Velocity  using  SFM  I  (RLM) 

The  velocity  model  corresponding  to  the  SFM  I  model  is  referred  to  as  the  refined  Langevin 
model  since  the  use  of  instantaneous  frequency  (rather  than  the  mean  frequency)  necessitates 
changes  to  the  simple  Langevin  model  (section  3.2)  as  described  in  Refs.  3-5  and  3-6.  The 
model  is  given  by 
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dU;  =- 


1  £dt  +g,dt-(|+|c0)(u*-<Ui>)(Y)<a»dt+didt+(Co0)*k')1'2dWi 


(3.6) 


The  salient  differences  between  the  simple  Langevin  model  SLM  I  (Eq.  3.3)  and  the  refined 
model  are  that  the  random  term  contains  the  instantaneous  frequency  (rather  than  the  mean),  the 
O)- weighted  mean  kinetic  energy  k'  is  used,  and  there  is  an  additional  drift  term  dj  which  is 
needed  to  account  for  the  correlations  between  U  and  0)  that  result  from  random  term.  The 
additional  drift  term  di  is  a  complicated  function  of  the  Reynolds-stress  tensor,  its  co-weighted 
counterpart  and  k',  as  discussed  in  Ref.  3-7.  Due  to  these  differences,  the  appropriate  value  of 
the  constant  for  C0  for  this  model  is  different  from  that  in  SLM  I.  The  recommended  value  is 
C0  =  3.5.  This  model  is  used  for  the  computations  reported  in  Chapters  7  and  8. 


3.5  Stochastic  Frequency  Model  II  (SFM II) 

The  stochastic  frequency  model  SFM  II  is  a  variation  of  the  SFM  I  model  and  is  the  most  recent 
model  developed  [3-10].  Compared  to  the  previous  stochastic  frequency  model  SLM  I  (section 
3.2),  the  new  model  is  easier  to  implement  and  is  expected  to  be  more  robust.  Not  only  is  the 
frequency  model  easier  to  implement  and  needs  fewer  computations,  but  also  the  corresponding 
velocity  model  is  much  simpler,  easier  to  implement  and  expected  to  be  more  robust. 

According  to  the  SFM  II  model,  the  increment  in  the  particle  frequency  over  an  infinitesimal 
time  interval  dt  is  given  by 

dca*  =-^-CcolSiiSiidt  -  m*<m>[Cco2+C3Qr]dt  +  <co>(2C3nrco*G2)1/2dW  +C3Q2dt  , 
<co>  J  J 

(3.7) 

where 

Qr  =  Q/<m>  ,  (3-9) 


and 


Q  =  Cn  <  0)  |  co*  >  <co>  > 


(3.10) 


3-4 


is  the  conditional  mean  frequency  defined  as  an  "above-average"  mean,  i.e.  it  is  proportional  to 
the  mean  of  the  instantaneous  frequencies  that  are  greater  than  or  equal  to  the  unconditional 
mean  frequency.  The  constant  Cq  (determined  in  terms  of  incomplete  gamma  functions)  has  the 
value  0.6893,  and  is  specified  such  that  Q  =<co>  in  homogeneous  turbulence. 

Although  equations  3.4  and  3.7  are  cast  in  the  same  general  form,  the  basis  for  the  derivation  of 
the  two  models  are  very  different,  and  the  corresponding  constants  appearing  in  the  equations 
have  significantly  different  values.  The  SFM  I  model  is  based  on  a  log-normal  distribution  for 
the  pdf  of  frequency  while  the  SFM  II  model  is  based  on  a  Gamma  distribution.  The  SFM  II 
model  was  also  constructed  with  reference  to  direct  numerical  simulation  results  [3-5].  By 
construction,  the  values  of  the  constants  C2  and  C3  are  0.25  and  1.0  respectively.  The 
recommended  values  for  and  Cm2  are  0.08  and  0.9  respectively. 

A  significant  feature  of  the  SFM  II  model  is  that  it  includes  the  conditional  mean  frequency 
which  will  be  significantly  higher  than  the  unconditional  mean  in  the  intermittent  region  and 
avoids  the  inclusion  of  an  ad  hoc  term  as  in  the  previous  model  to  account  for  intermittency  and 
entrainment  of  nonturbulent  fluid.  However,  a  drawback  of  the  model  compared  to  SFM  II  is 
that  the  conditional  (mean)  frequency,  and  not  the  instantaneous  particle  frequency,  appears  in 
the  particle  velocity  equation  corresponding  to  SFM  II  (see  next  subsection).  This  means  that  the 
full  effect  of  all  the  turbulence  time  scales  is  not  felt  by  the  velocity  equation. 

It  should  be  noted  that  both  the  SLM  I  and  SLM II  models  performed  equally  well  for  the  cases 
studied.  Even  though  SFM  II  may  be  simpler  to  implement  and  more  robust  than  SFM  I,  the 
latter  is  also  computationally  quite  tractable.  More  studies  will  be  required  before  a  clear 
preference  emerges. 


3.6  Simple  Langevin  Model  II  for  Velocity  using  SFM  II  (SLM  II) 

The  model  for  the  increment  in  the  particle  velocity,  that  corresponds  to  the  SFM  II  model,  is  a 
variant  of  the  simple  Langevin  model  (SLM  I)  described  in  section  3.2.  The  SLM  II  model  is 
described  by: 

jlj  1  /I  O  \  *  V 

dUi  ="7*aldt  +g'dt'  (2+4C°)(U|  -<ui>)<®>dt+(c.<®>k)^ 1  dWi  •  <3-3> 

As  explained  before,  the  first  two  terms  on  the  right  exactly  account  for  the  acceleration  due  to 
mean  pressure  gradients  (including  variable-density  effects)  and  body  forces  respectively.  The 
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last  two  terms  together  model  the  effects  of  viscous  dissipation  and  the  fluctuating  pressure 
gradient. 


The  main  difference  between  the  SLM  II  and  SLM  I  models  is  that  the  conditional  mean 
frequency  Q.  appears  instead  of  the  unconditional  mean  frequency  <co>.  The  consequences  of 
this  have  been  explained  in  the  previous  subsection.  Also,  the  value  of  C0  appropriate  for  the 
SLM  II  model  is  2.5.  The  SFM  II  and  the  SLM  II  models  have  been  used  to  compute  the  flows 
described  in  Chapter  10. 
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4.  MODELS  FOR  MOLECULAR  MIXING  OF  SCALARS 


This  chapter  describes  the  models  developed  and/or  used  in  the  present  program,  for  the 
molecular  mixing  of  scalars  in  the  particle  evolution  equations  (i.e.  the  term  involving  0  in  Eq. 
2-18).  Two  different  models  for  scalar  mixing  are  described  and  discussed.  One  of  the  models, 
the  improved  mixing  model,  is  a  "particle  interaction  model"  which  means  that  the  model  is 
based  on  the  direct  interaction  between  the  properties  (or  scalar  values)  of  the  individual  particles 
that  are  in  close  proximity.  The  other  is  a  relaxation  to  mean  model,  in  which  the  individual 
particles  interact  with  the  local  mean  value  of  the  scalar  which  is  of  course  based  directly  on  the 
properties  of  the  other  particles  in  the  vicinity. 


4.1  Improved  Mixing  Model  (IMM) 

The  improved  mixing  model  is  a  particle  interaction  model.  It  is  a  significantly  more 
sophisticated  and  improved  version  of  the  well  known  Curl’s  model  [4-1]  and  overcomes  many 
of  the  undesirable  features  of  the  Curl’s  model.  The  model  was  developed  by  Pope  [4-2]  and  has 
been  used  by  the  majority  of  joint  pdf  studies  as  well  as  for  nearly  all  the  flows  in  the  present 
study.  The  model  is  described  by  the  illustration  given  below. 

Consider  a  general  scalar  variable  <j).  In  the  mixing  model,  pairs  of  particles,  Np/2  in  number  (a 
total  of  Np  particles),  are  picked  are  picked  by  a  random  but  sophisticated  procedure  described  in 
Ref.  4-2.  The  number  Np  is  determined  by 

Np  =  q,  N  <co>  At  ,  (4.1) 

where  C^  is  a  model  constant  usually  in  the  range  1.5  to  2.0  (2.0  is  used  in  the  present  study),  N 
is  the  total  number  of  particles  in  the  computational  cell,  <co>  is  the  mean  turbulence  frequency 
in  the  cell,  and  At  is  the  interval  of  the  computational  time-step.  Let  p  and  q  denote  two  particles 
of  a  pair.  Their  scalar  values  change  from  those  at  time  t  to  those  at  time  t  +  At  as  follows: 

<|)p(t  +  At)  =  (!-£)  <|>p(t)  +  ,  (4.2) 


and 
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<j)q(t  +  At)  =  (1-5)  <|)q(t)  +  54>m  , 


(4.3) 


where 

♦m  =  H>P(t)  +  4>q(t)  ]  /  2  (4.4) 

is  the  common  mean  value  for  the  two  particles  and  ^  is  a  random  variable  with  the  distribution 
A(£)  (  =  10  [1  -  /4] ).  The  value  \  ( 0  <  \  <  1 )  controls  the  extent  of  the  mixing;  with  \ 

=  0  no  mixing  occurs,  with  \  =  1  maximum  mixing  occurs  and  both  particles  mix  to  their 
common  mean  4>m. 

Since  the  common  mean  <J)m  is  the  same  before  and  after  the  mixing  process,  the  local  mean 
scalar  value  <§>  calculated  as  an  ensemble  average  of  all  the  particles  in  the  cell  (section  2.2.2) 
does  not  change  due  to  mixing  process.  However,  since  the  mixing  process  brings  the  two 
mixing  particles  closer  to  each  other  in  the  (j)-space,  the  variance  of  the  .scalar  «()  2>  decreases. 
It  can  be  shown  that  the  effect  of  the  mixing  process  on  the  scalar  variance  is 

d«|>'  2>/dt  =  -  q,  <o)>  <$'  2>  •  (4-5) 


Thus,  the  mixing  model  serves  to  dissipate  the  fluctuations  in  the  scalar  and  is  used  to  model  the 
molecular  dissipation  process  represented  by  the  term  0  in  Eq.  2-18. 


4.2  Interaction  bv  Exchange  with  the  Mean  (EEM)  Model 

This  model  was  originally  suggested  by  Dopazo  [4-3].  It  is  a  deterministic  model  (i.e.  it  does  not 
involve  any  random  numbers),  and  is  only  based  on  the  scalar  value  of  the  particle  and  the  mean 
scalar  value  at  the  particle  location.  The  model  for  0  is  described  by; 

d<)>a*  =  -  Qj,  (<(>a*  -  «|>a>)  <co>  dt  .  (4.6) 


It  can  be  seen  that  the  model  has  the  desired  effect  of  moving  the  value  of  the  particle 
composition  towards  the  mean  composition  reducing  the  fluctuation  in  the  scalar  values.  Indeed. 
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the  variance  of  the  scalars  vary  as  in  Eq.  4.5.  The  value  of  the  constant  for  this  model  is  in  the 
range  1.0  to  1.5. 

The  model  is  very  simple  to  implement.  All  the  particles  mix,  and  there  is  no  need  to  develop  an 
algorithm  to  pick  the  pairs  of  paticles  to  mix.  Also,  for  general  variable-density  flows  additional 
care  should  be  taken  in  the  particle  interaction  models  to  ensure  that  the  mixing  is  weighted 
appropriately  to  account  for  the  differences  in  the  masses  of  the  mixing  particles.  But  with  the 
IEM  model,  the  mass  of  the  particle  is  not  a  consideration  in  carrying  out  the  mixing  process. 

However,  there  are  several  drawbacks  with  the  model.  Since  each  particle  changes  a  different 
amount  towards  the  mean  value  (depending  upon  its  own  scalar  value),  there  is  no  guarantee  that 
the  mean  value  «l>a>  remains  unchanged  after  the  mixing  process.  The  mean  is  conserved  in  a 
statistical  sense  only  when  there  are  a  large  number  of  particles  in  the  simulation.  Also,  if  the 
profiles  of  «j>a>  have  overshoots  and  undershoots  (over  and  below  their  allowed  values)  due  to 
the  nature  of  the  splines  used  for  smoothing  mean  fields,  the  particles  scalar  values  can  also  go 
out  of  bounds  as  a  result  of  this  mixing  model.  Another  major  deficiency  is  that  since  this  model 
considers  only  the  mean  value,  it  is  completely  missing  out  the  information  of  the  pdf  of 
composition  contained  in  the  ensemble  of  particles.  For  example,  since  it  does  not  use  the 
information  on  the  evolving  pdf,  starting  with  an  arbitrary  initial  pdf  of  scalars  in  a  mixing 
problem  in  homogeneous  turbulence,  this  model  will  not  relax  to  a  Gaussian  pdf  as  expected. 
Nevertheless,  due  to  its  extreme  simplicity  and  overall  good  performance  for  most  flows  with 
proper  initial  specification  of  the  pdf,  the  utility  of  this  model  for  complex  multi-dimensional 
flows  should  not  be  overlooked.  This  model  is  used  for  the  flows  computed  in  Chapter  10. 
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5.  ALGORITHMS  TO  DETERMINE  MEAN  PRESSURE 


This  chapter  briefly  describes  the  Monte  Carlo  solution  algorithm  and  the  algorithms  for 
determining  the  mean  pressure  within  the  solution.  Two  types  of  algorithms  are  described.  The 
first,  boundary-layer  algorithm,  is  applicable  to  one-dimensional  (1-D)  time-dependent  or  2-D 
planar  or  2-D  axisymmetric  flows  without  reverse-flows.  The  other  class  is  the  elliptic  flow 
algorithm  which  is  applicable  to  general  2-D  and  3-D  flows.  It  should  be  noted  that  the  flow  is 
1-D,  2-D  or  3-D  only  in  the  time-mean  sense.  In  all  the  pdf  calculations,  all  three  components  of 
physical  space  and  all  three  components  of  velocity  are  considered  in  the  computation  of  the 
particle  motions  and  velocities.  (Note:  For  the  discussion  presented  in  this  chapter  <P>  and  <p> 
refer  to  the  mean  pressure  and  mean  density,  respectively,  which  are  Reynolds-averaged 
quantities.  For  all  other  variables,  angled  brackets  <  >  indicate  Favre-averaged  quantities.) 


5.1  Boundary  Laver  Algorithm 

This  algorithm  is  applicable  to  parabolic  flows  (1-D  time-dependent  or  2-D  planar  or 
axisymmetric  flows  without  reverse-flow)  such  as  jets,  wakes,  mixing  layers,  etc.  The  boundary 
layer  assumptions  can  be  used  to  calculate  these  flows.  According  to  the  assumption  all  axial 
gradients  of  mean  quantities  are  assumed  to  be  small  compared  to  the  transverse  gradients  and 
are  neglected.  The  transverse  pressure  gradient  is  determined  from  the  solution  of  the  transverse 
mean  momentum  equation.  However,  in  the  present  study  involving  swirling  flows  (without 
recirculation)  the  axial  pressure  gradient  is  not  assumed  zero.  Instead  it  is  calculated  from  the 
transverse  (radial)  pressure  distribution  as  explained  later  in  this  section.  A  more  detailed 
explanation  of  the  material  in  this  section  can  be  found  in  Refs.  5-1  and  5-2. 

For  the  discussions  to  follow,  the  notation  Us(Uj,  U2,  U3)  =  (U,  V,  W)  =  (<U>  +  u,  <V>  +  v, 
<W>  +  w)  is  used.  In  addition,  due  to  the  geometry  of  the  flows  considered,  polar  cylindrical 
coordinates  (x,  r,  0)  are  used.  In  the  boundary-layer  algorithm,  the  solution  is  marched 
downstream  starting  from  initial  conditions.  The  joint  pdf  at  each  step  in  the  x-direction  (the 
predominant  flow  direction)  is  represented  by  a  large  number,  N,  of  computational  particles.  At 
the  axial  location  x,  each  particle  has  the  position  x*(x),  velocity  U*(x),  and  the  scalar  value 
<j)*(x).  In  the  general  step  from  x  to  x+Ax,  each  particle  evolves  over  at  time  interval  At*  given 
by: 


At*  =  Ax/Ui* 


(5.1) 
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Note  that  all  particles  advance  by  the  same  Ax,  but  have  different  time  intervals  At*.  The 
increment  to  their  velocities,  frequency  and  scalar  values  are  given  by  the  corresponding 
evolution  equations  discussed  in  the  previous  two  chapters. 


The  2-D  axisymmetric  case  applicable  to  swirling  flows  is  considered.  There  are  no  theoretical 
limitations  for  computing  swirling  flows  with  the  pdf  method.  Extensions  were  made  to  the 
boundary-layer  algorithm  to  be  able  to  compute  swirling  flows.  The  main  extensions  to  the 
algorithm  are  to  calculate  and  include  the  axial  pressure  gradient  and  to  include  all  terms  arising 
from  to  the  nonzero  <W>  field.  The  mean  radial  pressure  gradient  is  computed  from  the  mean 
radial  momentum  equation  with  the  boundary-layer  approximation  that  all  the  axial  gradients  are 
negligible.  The  azimuthal  gradients  are  zero  due  to  axisymmetry.  The  mean  radial  pressure 
gradient  is  given  by: 


1  3<P> 

<p>  3r 


1 

r 


^<W>2  +  <w  >2-<v>2-r 


3<v2>\ 
9r  )  ■ 


The  mean  pressure,  as  a  function  of  r  at  a  given  x,  is  given  by 


Ro 

fd<P> 

<P(r;x)>  =  <P(Rc;x)>  -  J — , 


(5.2) 


(5.3) 


where  Rq  is  the  radius  of  the  outer  boundary  of  the  jet.  The  mean  axial  pressure  gradient  at  x  + 
Ax,  as  a  function  of  r,  is  then  deduced  by 

3<P(r;  x+Ax)>  _  <P(r;  x+Ax)>-<P(r;  x)> 

3x  Ax 

Any  externally  imposed  axial  pressure  gradient  can  be  incorporated  through  <P(Ro;  x)>  in  Eq. 
5.3.  In  the  present  study,  <P(Rq;  x)>  is  taken  to  be  constant  for  all  x.  The  axial  pressure  gradient 
is  important  in  calculating  the  evolution  of  the  mean  axial  velocity  (especially  its  decay  on  the 
centerline)  when  appreciable  swirl  is  present  in  the  flow. 
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5.2  Mean  Pressure  Algorithm  I  for  Elliptic  Flows 


The  pressure  algorithm  is  constructed  for  steady-state,  constant-  or  variable-density,  high 
Reynolds  number  flows.  Details  of  the  algorithm  and  its  implementation  can  be  found  in  Ref. 
5-3.  The  equation  for  mean  pressure  is  the  Poisson  equation  (Eq.  5.6)  derived  from  the  mean 
momentum  equation  (Eq.  5.5): 


3<p>  <UjUj>  _  3<P> 

9xi  9x: 


(5.5) 


32<p>  32<p><uiuj> 


dxj  dx- 


0Xj  3xj 


(5.6) 


Since  the  right-hand  side  of  the  equation  can  be  evaluated  from  the  pdf  solution,  Eq.  5.6  can  be 
solved  for  the  mean  pressure  <P>  with  appropriate  boundary  conditions.  Starting  from  arbitrary 
initial  conditions,  if  an  iterative  or  psuedo-time  marching  algorithm  is  used  to  reach  the 
steady-state  for  variable-density  flows,  then,  because  the  transient  terms  are  absent  from  Eq.  5.6, 
the  mean  pressure  given  by  Eq.  5.6  is  not  correct  until  the  steady-state  is  achieved.  This  is  true 
even  for  the  constant  density  cases  unless  the  initial  velocity  field  satisfies  the  continuity 
equation  and  Eq.  5.6  is  solved  in  a  coupled  manner  with  the  mean  momentum  equation  which  is 
not  solved  directly  in  the  pdf  method.  A  consequence  is  that  the  mean  continuity  equation 


d<p>  <Uj> 
dxl 


=  0 


(5.7) 


will  not  be  satisfied.  Given  a  field  that  does  not  satisfy  Eq.  5.7,  a  velocity  correction  AU  can  be 
obtained  by 


AUt  = 


1  30 

<p>  3xj 


(5.8) 


By  requiring  that  {<p><U,>  +  <p>AU,}  be  divergence  free,  we  obtain  a  Poisson  equation  for  d>: 


52<j>  3<pxUi> 

3xj  3xj  3Xj 


(5.9) 
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These  observations  suggest  the  following  basic  algorithm:  every  so-many  steps 

1.  solve  Eq.  5.9  for  d> 

2.  add  the  velocity  correction,  Eq.  5.8 

3.  solve  Eq.  5.6  for  <P> . 

The  spirit  of  the  algorithm  is  similar  to  that  of  the  SIMPLER  algorithm  used  in  conventional  FV 
calculations  [5-4].  However,  the  solution  strategy  and  the  implementation  in  the  pdf  calculations 
are  markedly  different.  All  mean  quantities  in  the  pdf  method  are  calculated  by  forming 
ensemble  averages  and  smoothing  using  some  splining  procedure.  Since  the  Poisson  equation 
for  pressure  involves  second  derivatives  of  mean  quantities.  It  is  necessary  to  evaluate  the 
profiles  of  the  means  with  very  little  statistical  error.  Typically,  the  splines  used  are 
cross-validated  cubic  B-splines  [5-5]  in  order  to  obtain  the  required  accuracy.  But,  they  are  quite 
expensive  to  form  since  the  cross-validation  and  error  reduction  is  an  iterative  procedure. 

However,  an  efficient  and  economical  procedure  to  solve  the  set  of  equations  in  the  pressure 
algorithm  has  been  developed.  It  utilizes  the  spline  representation  of  the  mean  fields  and 
exploits  the  fact  that  the  resulting  matrix  equations  for  both  <P>  and  d>  have  the  same  form  and 
the  matrix  on  the  left-hand  side  of  the  two  equations  is  the  same  and  constant  through  out  the 
calculation.  Hence  the  expensive  matrix  inversion  is  performed  only  once  at  the  beginning  of  the 
calculations,  and  the  equations  for  <P>  and  O  are  solved  by  the  inexpensive  back-substitution 
procedure  on  subsequent  steps.  Additional  economies  can  also  be  realized  by  not  solving  these 
equations  on  every  step  but  once  in  a  few  steps  while  allowing  the  mean  velocity  fields  to 
change. 

Although,  this  algorithm  was  successfully  implemented  and  demonstrated  for  elliptic 
recirculating  flows  (see  Chapter  9),  it  is  felt  that  a  more  robust  and  economical  algorithm  will  be 
needed  for  more  complex  2-D  and  3-D  flows.  Also  it  would  be  very  cumbersome  and  difficult  to 
fit  cubic  B-splines  in  irregular  geometries  and  exploit  the  implementation  algorithm  described  in 
the  previous  paragraph. 


5.3  Mean  Pressure  Algorithm  II  for  Elliptic  Flows 

The  main  purpose  of  the  elliptic  flow  algorithm  is  to  determine  the  mean  pressure  field  to  be 
used  in  the  particle  velocity  equations  (Chapter  3)  while  ensuring  that  the  mean  conservation 
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equations  for  mass  and  momentum  are  satisfied.  The  elliptic  flow  algorithm  described  in  this 
section  is  an  alternative  to  the  one  described  in  the  previous  subsection.  In  as  much  as  the  role  of 
calculating  the  mean  pressure  is  to  ensure  that  the  mean  mass  and  momentum  equations  are 
satisfied,  this  algorithm  does  not  attempt  to  directly  solve  the  Poisson  equation  for  the  pressure 
but  performs  several  corrections  to  the  particle  properties  to  ensure  the  same  result.  The 
algorithm  performs  a  velocity  correction  to  satisfy  mean  mass  conservation  and  determines  a 
mean  pressure  correction  on  every  step  starting  from  arbitrary  initial  conditions.  Variance 
reduction  techniques  are  applied  so  that  mean  momentum  conservation  is  also  maintained  (i.e. 
turbulent  processes  such  as  mixing,  viscous  dissipation,  etc.  are  performed  on  subensembles  such 
that  the  subensemble  means  are  not  changed).  In  addition,  a  correction  to  the  position  of  the 
particles  is  made  to  ensure  that  the  consistency  condition  for  particle  methods,  namely  that  the 
volume  associated  with  a  subensemble  of  particles  should  equal  the  geometric  volume  occupied 
by  the  particles,  is  satisfied.  For  statistically  stationary  flows,  a  steady  state  is  achieved  in  which 
these  corrections  tend  to  zero  (in  the  mean,  and  the  variance  decreases  as  the  number  of  particles 
increases). 

In  the  algorithm  a  velocity  correction  potential  is  determined  such  that  after  adding  the 
velocity  correction 

AU  =  — —  V<D ,  (5.10) 

<P> 

(which  is  the  same  as  Eq.  5.8),  the  corrected  velocity  field  satisfies  the  continuity  equation  (Eq. 
5.7).  When  the  velocity  increment  is  determined  by  Eq.  5.10  for  a  time  step  At,  it  is  equivalent  to 
the  effect  of  a  mean  pressure  correction 

A<P>  =  -  <&/At  .  (5.11) 

The  Poisson  equation  for  the  velocity  correction  potential  (same  as  Eq.  5.9)  is  set  up  and  solved 
using  bilinear  basis  function  representation  for  calculating  mean  quantities  [5-7].  Thus,  the  mean 
pressure  field  is  not  determined  directly  from  the  solution  of  the  Poisson  equation.  However,  any 
error  in  the  mean  pressure  field  is  compensated  by  the  velocity  correction,  i.e.,  the  potential  O  is 
such  the  total  effect  of  the  correct  pressure  should  be  felt.  Computations  using  this  pressure 
algorithm  have  been  successfully  validated  against  data  for  nonreacting  and  reacting  flows  [5-8]. 
The  results  are  presented  in  Chapters  9  and  10. 

In  contrast,  the  pressure  algorithm  developed  and  used  by  Anand  et  al.  [5-3]  solves  for  the 
Poisson  equation  for  pressure  as  well  as  for  the  velocity  correction  potential.  However,  since  the 
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Poisson  equation  involves  second  derivatives  of  mean  velocities  it  is  necessary  to  determine  the 
mean  velocity  field  to  a  high  degree  of  accuracy.  Hence,  as  discussed  previously,  bidirectional 
cross-validated  cubic  splines  are  used  to  determine  means  in  that  algorithm  which  can  be 
computationally  expensive.  The  current  algorithm  is  expected  to  be  less  expensive  and  more 
robust.  The  more  important  advantage  is  that  it  is  easier  to  extend  the  current  algorithm  to 
irregular  geometries  (body-fitted  grids)  and  to  three-dimensional  flow  calculations. 
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6.  COAXIAL  JET  FLOWS 


6.1  Overview 


The  pdf  method  is  applied  to  calculate  the  flow  field  of  two  concentric  jets  issuing  into  stagnant 
surroundings.  The  main  objective  of  this  study  was  to  validate  the  pdf  method  using  the  mean 
frequency  model  (MFM)  and  the  simple  Langevin  model  (SLM  I)  for  general  inhomogeneous 
flows.  A  modeled  transport  equation  for  the  mean  turbulence  frequency  is  solved,  in  contrast  to 
previous  studies  with  the  pdf  method  that  assumed  the  frequency  to  be  constant  across  the  flow. 
This  results  in  a  significant  improvement  in  the  agreement  of  the  model  predictions  with 
experimental  data  of  McDonnel  et  al.  [6-1].  Other  aspects  of  the  modeling  are  investigated, 
including  the  modeling  of  intermittency,  the  effect  of  the  universal  constant  in  the  Langevin 
model  for  velocities,  and  the  assumption  of  a  coflow  for  the  simulation. .  The  study  is  reported  in 
detail  by  Anand  et  al.  [6-2] 


6.2  Background 

A  local  turbulent  time  scale,  i,  or  its  inverse,  the  turbulent  frequency,  <0)>,  is  needed  for  the 
models  use  in  the  pdf  method.  For  the  study  of  recirculating  flows  (Anand  et  al.  1987),  the 
turbulent  frequency  was  supplied  to  the  pdf  calculations  from  the  accompanying  conventional 
finite-volume  calculations.  In  all  other  studies  with  the  pdf  method,  the  turbulent  frequency  was 
assumed  to  be  constant  across  the  flow.  The  assumption  of  a  uniform  <co>  is  valid  by  definition 
for  the  homogeneous  turbulence  cases  [6-4,  6-5].  For  the  self-similar  free  shear  flows  [6-6,  6-7] 
the  turbulence  frequency  was  assumed  based  on  self-similarity  to  be  proportional  to  the  inverse 
of  the  mean-flow  time  scale.  The  same  assumption  was  used  for  calculating  the  <co>  in  the 
studies  of  jet  diffusion  flames  [6-8,  6-9]  where  the  emphasis  was  more  on  the  chemical  modeling 
aspects  than  on  the  velocity  fields.  The  assumption  of  a  uniform  <co>  is  perhaps  reasonable  in 
the  self-similar  regime  and  in  the  region  far  downstream  of  the  inlet  for  the  free  shear  flows 
studied  previously.  However,  the  effect  of  a  uniform  <co>  on  the  pdf  calculation  in  the 
developing  (non-self-similar)  region  close  to  the  inlet  or  for  other  inhomogeneous  flows  need  to 
be  investigated. 
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6.3  Modeling  Considerations 


The  boundary  layer  algorithm  (section  5.1)  is  used  since  the  flow  is  parabolic  in  nature.  The  mea 
frequency  model  (MFM)  (section  3.1)  is  used  for  the  frequency  and  simple  Langevin  model 
(SLM  I)  is  used  for  the  velocities.  Calculations  were  also  made  with  the  assumption  of  uniform 
<G»  for  comparison. 

In  the  boundary-layer  algorithm  (see  section  5.1),  modeled  particles  representing  fluid  particles 
are  marched  downstream  based  on  their  axial  velocities.  It  is  therefore  necessary  that  the  axial 
velocities  be  always  positive.  Hence,  a  coflow  stream  with  a  small  positive  axial  velocity 
surrounding  the  jets  is  assumed  for  the  simulations,  although  the  jets  issue  into  stagnant 
surroundings  in  the  experiments.  The  effect  of  this  assumed  coflow  is  also  investigated. 

Other  aspects  of  the  modeling  that  are  investigated  in  the  present  study  are  the  modeling  of 
intermittency  and  the  effect  of  the  universal  constant  in  the  Langevin  model  for  the  velocities. 

6.3.1  Intermittency 

The  phenomenon  of  intermittency,  which  results  from  nonturbulent  fluid  surrounding  the 
turbulent  flow  under  consideration,  is  accounted  for  by  conditional  modeling  based  on  a  passive 
scalar  [6-6,  6-8].  The  value  ())*>0  represents  turbulent  fluid  and  (J)*=0  denotes  nonturbulent 
fluid.  Within  the  turbulent  fluid,  the  (j>*  value  of  a  particle  evolves  by  molecular  mixing,  which 
is  modeled  by  the  IMM  model  (section  4.1). 

In  the  intermittency  model,  nonturbulent  fluid  is  entrained  (i.e.,  becomes  turbulent)  at  a  specific 
rate.  The  nonturbulent  fluid  exchanges  momentum  and  energy  with  the  turbulent  fluid.  The 
model  involves  three  constants:  Cg  (=2.0)  controls  the  rate  of  entrainment;  Cm  (=1.5)  and  Ce 
(=5.0)  control  the  rate  of  momentum  and  energy  transfer,  respectively.  Haworth  and  Pope  [6-6] 
concluded  that  although  the  intermittency  model  works  reasonably  well  in  conjunction  with 
particle  interaction  models  for  velocity  (not  used  here),  it  is  not  entirely  satisfactory  in 
conjunction  with  the  Langevin  model  which  is  the  preferred  model.  In  the  present  study, 
calculations  are  performed  both  with  and  without  the  intermittency  model. 

6.3.2  Calculations  with  Uniform  <co> 

For  the  calculations  in  the  present  study  for  which  <to>  is  assumed  constant  in  the  r-direction,  the 
last  term  (diffusion  term)  in  Eq.  3.1  is  zero  by  assumption.  The  rest  of  the  terms  in  the  equation 


6-2 


are  weighted  by  the  local  turbulent  kinetic  energy  and  integrated  over  the  cross-sectional  area  of 
the  flow  to  obtain  an  integral  estimate  for  the  value  of  <co>. 


6.4  Flow  Description  and  Computational  Details 

The  flow  studied  is  of  two  concentric  nonswirling  jets  issuing  into  stagnant  surroundings. 
Detailed  experimental  data  for  the  flow,  including  mean  velocities  and  turbulent  normal  and 
shear  stresses,  are  reported  in  Ref.  6-1.  The  central  air  jet  of  diameter  D  (24.1  mm)  is 
surrounded  by  an  annular  air  jet  with  inner  and  outer  diameters  of  29  mm  and  36.7  mm, 
respectively.  The  average  axial  velocity  in  the  central  jet  is  4.67  m/s  and  that  in  the  annular  jet  is 
6.6  m/s.  The  Reynolds  number  of  the  central  jet  (based  on  D)  is  approximately  6100. 

Initial  conditions  for  the  computations  are  prescribed  from  experimental  data.  Due  to  the  wall 
thickness  of  2.5  mm  between  the  central  and  annular  jets,  a  reverse  flow  region  is  indicated  by 
the  data  at  x/D=0.08  form  the  nozzle  exit.  Hence,  the  next  axial  measurement  station  at 
x/D=0.62  is  taken  as  the  initial  plane  for  the  calculations. 

In  this  and  the  next  section,  the  following  notation  is  used:  <U>  is  the  mean  axial  velocity;  u,  v, 
and  w  are  the  fluctuating  velocities  in  the  axial,  radial,  and  azimuthal  directions.  The  reference 
velocity  used  for  normalization  (Uox  =  4.28  m/s)  is  the  mean  axial  velocity  on  the  centerline  at 
the  initial  station.  The  radial  distance,  r,  is  normalized  by  the  central  jet  radius,  R  (=D/2). 

The  initial  velocity  pdf  is  prescribed  to  be  joint  normal  with  the  mean  and  covariances  taken 
from  linearly  interpolated  experimental  data  available  (see  Figure  6-1).  Two  sets  of 
experimental  data  available  for  <U>  and  <u2>  give  an  indication  of  the  scatter  in  the  data.  The 
nominal  number  of  particles  used  for  the  computations  is  23,000.  A  coflowing  stream  of  mean 
axial  velocity  Ub=0.4  m/s  is  assumed.  With  the  initial  peak  mean  axial  velocity  being 
approximately  6  m/s,  the  velocity  ratio  (low  speed  to  high  speed)  is  0.066. 


The  initial  values  of  <co>  are  derived  from  experimental  data  based  on 


<(0> 


,  k  3<U> 
^  <uv>  dr 


(6.1) 


which  follows  from  the  k-e  model  assumptions.  The  experimental  data  for  k,  <uv>,  and  <U>  are 
smoothed  using  B-splines  to  evaluate  Eq.  6.1. 
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The  initial  profile  of  <co>  from  Eq.  6.1,  the  integrated  uniform  value  estimated,  and  the  initial 
profile  of  <(j)>  are  shown  in  Figure  6-2. 

The  initial  pdf  of  <|)  was  taken  to  be  a  double-delta  function  of  magnitude  «[»  at  \[r=1.0  and 
magnitude  (1-«|)>)  at  \|/=0  .  Hence,  initially,  the  profile  of  «]»  corresponds  to  the  profile  of  the 
intermittency  factor  y,  which  is  defined  as  the  probability  of  the  fluid  being  turbulent.  The 
profile  of  «(»  is  specified  to  vary  smoothly  from  <4»=1  to  «j»=0  in  the  region  where  the 
kinetic  energy  profile  drops  in  the  outer  shear  region  of  the  annular  jet. 
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Figure  6-1.  Initial  conditions  for  normalized  mean  axial  velocity  and  Reynolds  stresses. 
Symbols  represent  experimental  data;  solid  lines  are  calculated  from  initial  pdf.  <w2>  is  similar 

to  <v2>. 


6-4 


Figure  6-2.  Initial  profiles  of  turbulence  frequency  <co>  and  mean  of  the  conserved  scalar  «|». 


6.5  Results  and  Discussion 


Limited  results  are  presented  first  for  three  cases  with  variations  in  the  modeling  used.  Case  A 
includes  modeling  of  intermittency  and  the  assumption  of  a  uniform  <co>.  For  Case  B, 
intermittency  is  not  modeled  (i.e.,  the  entire  flow  is  assumed  turbulent),  but  <oo>  is  assumed  to 
be  uniform.  For  Case  C,  intermittency  is  not  modeled,  and  <co>  is  nonuniform  (MFM  model). 

Results  for  the  mean  axial  velocity,  turbulent  kinetic  energy,  and  shear  stress  for  the  three  cases 
are  shown  in  Figure  6-3  at  two  axial  stations.  Overall,  Case  C  has  the  best  agreement  with  data, 
while  Case  A  has  the  worst.  The  difference  between  modeling  (A)  and  not  modeling 
intermittency  (B)  affects  the  peak  and  outer-edge  values  of  turbulent  stresses.  However,  the  jet 
spreading  rate  for  both  cases  is  smaller  compared  with  data,  although  not  modeling  intermittency 

(B)  improves  the  spreading  rate  and  peak  values  of  turbulent  stresses. 

A  more  significant  improvement  in  the  results  is  seen  when  <to>  is  assumed  to  be  nonuniform 

(C) .  The  mean  axial  velocity  profile  spreads  more  in  accordance  with  the  data,  while  the 
magnitude  of  turbulent  kinetic  energy  has  increased  in  the  central  and  outer  regions  of  the  jet. 
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Figure  6-3.  Comparison  of  results  for  cases  (A)  with  conditional  modeling  and  uniform  <co>, 
(B)  without  conditional  modeling  and  uniform  <CD>,  and  (C)  without  conditional  modeling  and 
nonuniform  <co>,  at  two  axial  stations.  Symbols  represent  experimental  data. 


This  is  because  a  higher  value  of  <co>  for  Cases  A  and  B  than  for  Case  C  in  these  regions 
(Figures  6-2  and  6-6)  implies  a  higher  ratio  of  dissipation  to  kinetic  energy  for  the  former  cases. 

Based  on  these  results,  the  modeling  strategy  adopted  for  the  rest  of  the  study  is  to  solve  the  full 
transport  equation  for  <co>  and  not  to  use  conditional  modeling  as  in  Case  C.  Detailed  results  for 
Case  C  and  comparison  with  data  are  presented  in  Figures  6-4  and  6-5  at  various  axial  stations  up 
to  x/D=12.4  (Results  and  data  for  <w2>  are  qualitatively  and  quantitatively  similar  to  those  for 
<v2>.)  The  computed  results  are  generally  in  good  agreement  with  data.  The  calculated  mean 
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Figure  6-4.  Computed  (Case  C)  radial  profiles  of  normalized  mean  axial  velocity  and  turbulent 
kinetic  energy  compared  against  experimental  data  at  various  axial  stations. 
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Figure  6-5.  Computed  (Case  C)  radial  profiles  of  normalized  Reynolds  stresses  compared 
against  experimental  data  at  various  axial  stations. 


axial  velocity  profiles  and  the  spreading  of  the  jet  are  well  predicted,  although  the  mean  axial 
velocity  in  the  central  region  of  the  jet  decays  slightly  faster  than  indicated  by  the  data.  The 
shear  stress  profiles  are  in  good  agreement  with  the  data.  Results  from  normal  stresses  are 
discussed  later  in  this  section. 

The  calculated  profiles  of  <co>  for  Case  C  are  shown  in  Figure  6-6.  The  peak  value  of  <to> 
drops  dramatically  from  about  700  s'1  at  a  location  less  than  half  a  diameter  downstream 
(x/D=1.04,  Figure  6-6).  The  <co>  profile  becomes  progressively  flatter  and  lower  in  magnitude 
further  downstream. 

The  effect  of  the  value  of  in  the  MFM  model  (Eq.  3.1)  is  now  investigated.  Also  shown  in 
Figure  6-6  are  <to>  profiles  with  affl=0.7  and  aw=1.3,  as  well  as  the  values  of  uniform  <co> 
calculated  in  Case  B.  Over  most  of  the  radius,  the  calculated  uniform  values  of  <G)>  are  much 
higher  than  those  calculated  in  Case  C.  It  is  also  apparent  that  the  value  of  has  little  effect  on 
the  results.  This  indicates  that  the  diffusion  term  in  the  <co>-equation  (Eq.  7)  is  small  compared 
to  the  other  terms  in  the  equation  for  the  flow  being  studied.  The  improvement  in  the  results  for 
Case  C  over  Case  B  comes  not  from  just  including  the  diffusion  terms,  but  more  due  to  the  fact 
that  a  local  value  of  <co>  is  calculated  rather  than  using  a  uniform  integral  value. 


co/100  ($-1) 


Figure  6-6.  Computed  radial  profiles  of  <co>  for  Case  C  (a^  =  1.0),  Case  B  (uniform  <co>),  and 
the  two  other  values  of  aw  compared  at  different  axial  stations. 
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A  discussion  of  the  calculated  kinetic  energy  and  velocity  variances  is  now  presented.  Figures 
6-4  and  6-5  show  that  while  the  k  profiles  are  in  good  agreement  with  the  data,  <u2>  is 
overestimated,  while  <v2>  (and  <w2>)  are  underestimated,  especially  at  the  near-inlet  stations. 
Since  the  redistribution  of  k  is  influenced  by  C0,  calculations  were  performed  with  different 
values  of  C0  in  the  range  2.1  to  3.5.  A  larger  value  of  C0  is  expected  to  increase,  as  one  of  its 
effects,  the  tendency  towards  isotropy  of  normal  stresses.  Results  show  that  as  C0  is  increased 
up  to  3.0,  the  agreement  of  <u2>,  <v2>,  and  <w2>  with  data  does  not  improve  up  to  x/D=3.11. 
At  x/D=6. 1 1  and  beyond,  the  results  are  independent  of  C0.  Also  the  results  at  all  x/D  are  nearly 
identical  for  Co=3.0  and  C0=3.5.  Overall,  the  effect  of  varying  C0  is  not  significant  enough  to 
warrant  a  revision  of  the  value  C0=2.1  used  in  this  and  previous  studies.  It  should  be  mentioned 
that  the  more  complicated  generalized  Langevin  model  [6-5],  which  explicitly  incorporates 
Reynolds  stresses  and  mean  velocity  gradients  in  the  model,  may  be  more  accurate  in  calculating 
the  relative  magnitudes  of  the  normal  stresses  than  the  simplified  model  used  here. 

Lastly,  the  effect  of  the  assumed  coflow  is  investigated.  Figure  6-7  shows  the  mean  axial 
velocity  profiles  for  three  cases  with  values  of  Ub  equal  to  0.4,  0.6,  and  0.8  m/s.  It  is  seen  that 
the  velocity  profiles  are  nearly  identical  except  near  the  edge  of  the  flow.  The  value  of  Ub  had 
no  significant  effects  on  the  turbulent  stresses.  Thus,  over  the  domain  of  interest  in  the  present 
flow,  the  value  of  Ub  has  little  effect  on  the  computed  results. 

The  improvement  in  the  results  seen  with  nonuniform  <0)>  in  the  present  study  can  be  viewed 
from  another  perspective.  With  similar  modeling  but  uniform  <co>,  Haworth  and  Pope  [6-7] 
found  that  the  boundary-layer  algorithm  was  increasingly  inaccurate  and  the  spreading  rates  were 
significantly  underestimated  for  values  of  the  low-speed  to  the  high-speed  velocity  ratio  less  than 
0.3.  The  results  in  the  present  study  with  a  velocity  ratio  of  0.066  indicate  that  the  range  of 
applicability,  in  terms  of  the  velocity  ratio,  is  extended  for  the  boundary-layer  algorithm  when 
the  local  value  of  <co>  is  used  instead  of  an  integral  value. 

An  additional  advantage  derived  in  solving  for  nonuniform  <co>  has  to  do  with  the  computational 
time  required.  The  distance  that  the  solution  can  be  marched  on  each  step  (Ax)  is  determined 
*  primarily  by  the  local  turbulent  time  scale.  The  average  At*  for  particles  is  limited  to  a  fraction 
of  the  local  value  of  x  (or  l/<co>).  From  Eq.  1,  it  can  then  be  seen  that  Ax  is  proportional  to  the 
minimum  value  of  <U>/<co>.  For  the  cases  in  the  present  study  with  nonuniform  <co>,  regions 
of  low  <U>  have  low  values  of  <co>  also,  while  for  the  uniform  cases. 


6-10 


«J>/U0,c 


Figure  6-7.  Mean  axial  velocity  profiles  from  calculations  with  three  different  values  of  coflow 

velocity  U5. 


the  value  of  <co>  is  significantly  higher  in  these  regions.  Hence,  the  allowable  Ax  for  each  step 
is  lower  for  the  uniform-<co>  cases  than  for  the  nonuniform  cases.  The  computations  were 
performed  on  an  AMDAHL  5890  mainframe.  For  the  nonuniform  <co>  cases,  330  steps  (5.8  min 
of  CPU  time)  are  typically  needed  for  the  calculations  from  x/D=0.62  to  x/D=12.4.  The  cases 
with  uniform  <co>  required  nearly  60%  more  steps  and  CPU  time. 


6.6  Summarizing  Remarks 

A  detailed  modeling  study  of  coaxial  jets  issuing  into  stagnant  surroundings  is  reported.  A 
modeled  transport  equation  for  the  turbulence  frequency  (<co>)  is  proposed  and  solved  to  obtain 
the  local  turbulent  time  scale.  The  results  of  these  calculations  are  in  better  agreement  with  the 
data  than  with  the  assumption  of  a  uniform  <co>  across  the  flow.  The  former  calculations  also 
require  significantly  less  computational  time  than  the  latter  to  advance  the  solution  by  the  same 
axial  distance.  The  simplified  Langevin  model  with  the  value  of  the  Kolmogorov  constant 
C0=2.1  was  found  to  be  satisfactory,  although  increasing  C0  produced  better  results  in  a  limited 


sense.  The  modeling  of  intermittency  was  found  to  be  deficient.  The  assumed  coflow  needed 
for  the  simulations  had  little  effect  on  the  results. 

Some  of  the  modeling  issues  raised  in  the  present  study  will  be  addressed  by  the  new  velocity- 
frequency  models  described  in  sections  3-3  through  3-6.  With  these  models,  the  turbulent 
frequency  is  also  considered  a  random  variable  similar  to  the  velocities  and  is  solved  as  a  part  of 
the  pdf  solution.  Calculations  using  these  models  are  described  in  the  following  chapters. 
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7.  CONSTANT-DENSITY  SWTRLING  JET  FLOWS 


7.1  Overview 

There  are  three  aspects  to  the  study  presented  in  this  chapter:  i)  swirling  flows,  an  essential 
ingredient  in  practical  combustion  devices,  had  not  been  studied  previously  with  the  pdf  method. 
In  the  present  study,  the  pdf  method  is  extended  and  applied  to  a  coaxial  swirling  jet  flow  and  the 
results  are  compared  with  detailed  experimental  data;  ii)  the  newly  developed  stochastic 
frequency  model  (SFM  I,  section  3.3),  initially  developed  and  validated  for  homogeneous 
turbulence  and  self-similar  flows  [7-1,  7-2],  is  applied  to  and  validated  in  conjunction  with  the 
refined  Langevin  model  (RLM)  for  a  general  developing  shear  flow  as  represented  by  the  coaxial 
swirling  jet  flow  considered.  The  results  are  also  compared  against  computations  using  the 
MFM  and  SLM I  models  (see  Chapter  3);  iii)  detailed  velocity  statistics,  conditional  upon  the  jet 
from  which  the  fluid  originated,  are  calculated  and  compared  .against  measurements, 
demonstrating  that,  in  addition  to  other  advantages,  the  pdf  method  can  provide  considerably 
more  description  of  the  flow  than  conventional  models,  and  yet  is  computationally  tractable.  The 
results  are  in  very  good  agreement  with  data.  The  experiments  were  jointly  designed  by  Allison, 
UDRI  and  WPAFB  personnel.  The  experiments  were  set  up  and  conducted  by  UDRI  personnel 
at  WPAFB.  A  detailed  account  of  the  study  presented  in  this  chapter  can  be  found  in  Ref.  7-2. 


7.2  Background 

Swirling  flows  (flows  with  appreciable  bulk  swirl)  are  an  essential  feature  in  many  practical 
combustors  such  as  gas  turbine  combustors.  Swirling  flows  had  not  been  studied  previously  with 
the  pdf  method,  although  there  are  no  theoretical  restrictions  for  applying  the  pdf  method  for 
such  flows.  For  elliptic  flow  calculations,  the  treatment  of  swirl  is  automatically  included  as 
long  as  no  assumption  is  made  that  the  mean  swirl  velocity  is  zero  and  all  the  terms  involving  the 
mean  swirl  velocity  are  retained.  In  the  present  study,  the  necessary  extensions  to  the  boundary 
layer  solution  algorithm  are  made  to  include  the  treatment  of  swirl  (see  section  5.1).  A  salient 
feature  of  this  extension  is  the  inclusion  of  the  mean  axial  pressure  gradient  in  the  calculations. 

It  is  well  known  that  the  k-e  model,  which  is  the  most  widely  used  turbulence  model,  is 
unsatisfactory  for  swirling  flows  due  to  effects  of  streamline  curvature,  although  the  model 
performs  reasonably  well  for  nonswirling  jet  flows.  Several  modifications  to  the  k-e  model  to 
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account  for  curvature  effects,  through  corrections  based  on  the  Richardson  number,  have  been 
investigated  and  reported  in  the  literature.  The  corrections  are  typically  made  through  the 
equation  for  mean  dissipation  as  a  means  of  correcting  the  turbulent  transport.  However,  these 
modifications  have  not  proved  to  be  uniformly  successful  over  a  range  of  swirling  flows.  In  this 
context,  it  is  necessary  to  verify  the  performance  of  the  pdf  method  for  swirling  flows. 

The  stochastic  frequency  model  is  used  for  the  computations.  This  represents  the  one  of  the  first 
applications  of  the  joint-velocity-frequency-scalar  approach  for  a  general  inhomogeneous  flow. 
In  addition  to  effecting  complete  closure  of  the  joint  pdf  equation,  the  inclusion  of  the 
instantaneous  frequency  in  the  joint  pdf  leads  to  more  realistic  modeling.  For  example,  as 
discussed  in  section  3.3,  the  model  allows  for  multiple  time  scales  rather  than  a  single  mean  time 
scale,  accounts  for  internal  intermittency,  has  potential  to  account  for  large  scale  structures,  and 
better  accounts  for  the  influence  of  the  origin  and  history  of  the  fluid  particles  on  the  local 
turbulent  structure. 

The  pdf  calculations  are  validated  against  benchmark  quality  data  of  Takahashi  et  al.  [7-4]  for  a 
coaxial  swirling  jet  flow.  A  unique  feature  of  the  laser  Doppler  velocimeter  (LDV) 
measurements  of  Takahashi  et  al.  [7-4]  is  that  conditional  velocity  statistics  are  measured  by 
seeding  individual  jets  one  at  a  time.  This  eliminates  issues  of  velocity  bias  in  the 
measurements.  Conditional  turbulent  correlations  up  to  fourth  order  have  been  reported.  These 
conditional  correlations  are  easily  computed  in  the  pdf  method,  and  the  comparisons  against  data 
serve  to  verify  that  turbulent  transport  is  indeed  accurately  calculated  in  the  pdf  method. 


7.3  The  Joint  PDF  Method  and  Models  Used 

The  boundary-layer  algorithm  described  in  section  5.1  with  the  extension  for  swirling  flows  is 
used.  Both  the  mean  and  stochastic  frequency  models  (MFM  and  SFM  I)  are  used.  The 
corresponding  velocity  models  (SLM  I  and  RLM)  are  used.  (See  Chapter  3  for  details  of  the 
models.)  The  MFM  model  is  also  referred  to  as  the  <co>-model  in  discussions  in  the  present 
chapter. 
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7.4  Flow  Description  and  Computational  Details 


7.4.1  Flow  Description 

The  flow  considered  has  a  central  nonswirling  jet  of  diameter  D  =  9.45  mm,  a  concentric  annular 
swirling  jet  of  diameter  Da  =  26.92  mm  surrounded  by  a  low  velocity  coflow  in  a  test  section  of 
150  mm  diameter.  A  schematic  of  the  swirling  jet  diffusion  flames  set  up  and  details  of  the 
diagnostics  used  are  presented  in  the  next  chapter.  For  the  experiments  used  in  this  chapter,  all 
the  jets  are  of  air.  The  bulk-averaged  axial  velocities  for  the  central  jet,  annulus,  and  coflow  are 
100,  20,  and  4  m/s,  respectively.  The  Reynolds  number  for  the  central  jet  is  approximately 
21,400. 

The  swirl  in  the  annular  jet  is  produced  by  a  vaned  swirler  with  a  vane  helix  angle  of  30  deg, 
located  96  mm  upstream  of  the  nozzle  exit.  The  swirl  number  for  the  annular  jet,  Sa,  calculated 
from  measured  data  1.5  mm  downstream  of  the  nozzles,  is  0.41.  However,  the  overall  swirl 
number  Sov,  based  on  the  central  and  annular  jets  (not  including  the  coflow)  is  0.09  due  to  the 
large  bulk  axial  velocity  of  the  central  jet.  Since  the  overall  swirl  number  is  low,  the  inclusion  of 
the  axial  pressure  gradient  in  the  boundary-layer  algorithm  did  not  have  a  significant  effect  on 
the  results.  However,  with  another  test  problem  [7-5]  with  an  overall  swirl  number  of 
approximately  0.4,  the  inclusion  of  the  axial  pressure  gradient  was  necessary  to  obtain  the  correct 
(measured)  spreading  of  the  mean  axial  velocity  profile.  The  latter  results  are  not  shown  here, 
but  are  reported  in  Ref.  7-6. 


7.4.2  Initial  Conditions 

Initial  conditions  for  the  computations  are  prescribed  from  experimental  data.  The  first 
measurement  station  x  =  1.5  mm  (x/D  =  0.16)  is  taken  as  the  initial  plane  for  the  calculations. 
The  initial  velocity  pdf  is  prescribed  to  be  joint  normal  with  the  mean  and  covariances  taken 
from  linearly  interpolated  experimental  data.  The  initial  pdf  of  frequency  for  the  stochastic 
model  is  taken  to  be  log-normal,  in  accordance  with  the  construction  of  the  model.  The 
specification  of  the  initial  <to>  profile  for  the  stochastic  model  is  discussed  later  in  this  section. 

The  initial  profile  of  <co>  for  the  <co>-model  is  derived  from  experimental  data  using  the 
expression  based  on  the  k-e  model: 
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(7.1) 


where  the  equivalent  velocity  gradient  and  shear  stress  are  given  by: 


and 

<uv>eq  =  (<uv>2  +  <uw>2) 1/2 

The  initial  <co>  profile  was  also  specified  from  Eq.  7.1  with  the  assumption  of  local  equilibrium 
of  turbulence  energy  in  the  shear  layer  that  yields  <uv>eq  =  C^2  k  so  that: 


(7.2) 

(7.3) 


<CD>  = 


(7.4) 


Although  Eqs.  7.1  and  7.4  result  in  slightly  differing  profiles  of  initial  <co>,  both  specifications 
produced  nearly  the  same  results  at  the  downstream  stations.  For  the  results  with  the  <co>-model 
reported  here,  the  specification  in  Eq.  7.4  was  used. 

The  specification  of  <co>  for  the  stochastic  model  needs  some  clarification.  The  production  of 
<co>  in  the  stochastic  model,  controlled  by  Cwi  (=  0.04)  (see  Eq.  3.3),  is  only  40%  of  the 
production  with  the  <to>  model  [Eq.  3.1,  Cj  =  0.1].  The  dissipation  of  <co>  is  also  lower  in  the 
stochastic  model,  but  is  about  80%  of  that  with  the  <co>-model  [C^  =  0-70  versus  C2  =  0.88]. 
Hence,  with  the  two  models  producing  nearly  the  same  velocity  statistics,  the  stochastic  model 
will  yield  lower  values  of  <to>,  especially  in  the  near-nozzle  region  where  production  is 
dominant  due  to  large  velocity  gradients.  This  fact  is  borne  out  by  the  results  presented  in  the 
next  section.  With  the  specification  of  the  same  inlet  <to>  for  both  the  models,  the  stochastic 
model  yielded  lower  <co>  (by  about  30%)  than  the  <co>-  model  at  the  downstream  locations  (x/D 
=  0.5  to  10  approximately).  The  agreement  with  data  of  the  turbulent  correlations  close  to  the 
nozzle  (x/D  =  1.06  and  2.65  stations)  was  poorer.  In  particular  the  turbulent  kinetic  energy  was 
lower  than  in  the  experiments  in  this  region.  This  would  have  an  effect  on  the  results 
downstream,  especially  the  spreading  rate.  The  results  suggested  that  the  appropriate  value  of 
inlet  <co>  for  the  stochastic  model  is  about  70%  of  the  value  used  for  the  <co>  model,  keeping  in 
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mind  that  the  stochastic  model  has  fluctuations  in  co  and  noting  the  differences  in  the  velocity 
models  discussed  previously.  Hence  the  inlet  <co>  for  the  stochastic  model  is  specified  to  be  0.7 
times  the  value  obtained  from  Eq.  7.4.  Again,  the  alternative  specification  of  0.7  times  the  value 
obtained  from  Eq.  7.4  produced  nearly  the  same  results. 

To  utilize  the  opportunity  provided  by  the  stochastic  model  to  specify  additional  initial 
conditions  based  on  joint  velocity-relaxation  rate  statistics,  initial  conditions  on  <v  co>  were 
specified  such  that  the  correlation  coefficient  for  <v  co>  is  the  negative  (from  wall 
boundary-layer  considerations)  of  that  of  <uv>.  Although,  this  specification  did  not  have  a 
major  effect  on  the  results,  it  underscores  the  fact  that  more  information  is  contained  in  and  can 
be  supplied  to  the  stochastic  model. 


7.4.3  Conditional  Statistics 

As  mentioned  previously,  LDV  measurements  for  the  flow  were  made  by  seeding  individual  jets 
one  at  a  time.  The  complete  flow  field  would  be  mapped  by  seeding  a  particular  jet  and  the 
procedure  would  be  repeated  for  the  remaining  jets.  The  experiments  were  shown  to  be 
repeatable.  This  technique  measures  the  conditional  pdfs  (i.e.,  conditional  upon  the  jet  of  origin) 
and  avoids  errors  due  to  velocity  bias  that  will  be  significant  due  to  the  large  differences  in  the 
bulk  velocity  among  the  jets.  However,  the  unconditional  pdf  cannot  be  deduced  from  these 
measurements  without  the  additional  measurement  of  the  fraction  of  time  the  fluid  from  each  jet 
spends  within  the  sample  volume.  The  advantage  with  the  pdf  method  is  that  the  conditional 
pdfs  are  easily  calculated  in  addition  to  the  unconditional  pdf.  Conventional  turbulence  models 
can  only  calculate  (a  limited  number  of)  unconditional  correlations  and  cannot  calculate 
conditional  quantities  without  requiring  additional  modeling. 

The  comparison  between  conditional  measurements  and  calculations  sheds  valuable  light  on  the 
detailed  performance  on  the  model  and  on  the  dynamics  of  individual  jets.  In  the  present  study, 
good  agreement  between  calculated  and  measured  conditional  quantities  also  implies  that  the 
calculated  unconditional  quantities  are  correct  since  no  additional  modeling  is  required  to  derive 
the  conditional  quantities. 

In  the  calculations,  the  particles  are  tagged  at  the  initial  plane  by  assigning  the  value  of  c*  to  be 
1,  2,  or  3  depending  on  their  jet  of  origin.  The  value  of  c*  does  not  change  throughout  the 
calculation.  Conditional  quantities  are  computed  for  output  by  forming  sums  over  particles  with 
the  same  value  of  c*  within  a  spatial  cell  and  fitting  B-splines  to  these  sums. 
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7.4.4  Computational  Parameters 


The  total  number  of  particles  used  in  the  simulations  is  150,000.  A  large  number  is  used  to 
reduce  statistical  error  especially  in  the  calculation  of  the  third  and  fourth  order  moments  (mean 
and  second  order  quantities  can  be  calculated  to  within  a  few  percent  accuracy  with  as  few  as 
25,000  particles).  The  number  of  cells  used  for  forming  sums  is  44  and  the  number  of  basis 
functions  used  for  B-splines  is  20. 

The  stochastic  model  calculations  took  approximately  38  minutes  of  CPU  time  to  march  from 
x/D  =  1.06  to  x/D  =  26.5  (in  about  540  steps)  on  a  CRAY  YMP.  The  <G»-model  took  about  28 
minutes  and  470  steps.  No  significant  effort  has  been  made  to  optimize  the  computer  code  for 
speed.  The  storage  required  was  about  7  megawords. 


7.5  Results  and  Discussion 

The  measured  mean  axial  velocity  on  the  centerline  of  the  jet  at  the  x  =  1.5  mm  station,  <U>0C 
(=129.7  m/s),  is  used  for  normalizing  all  the  velocity  statistics  presented.  The  radial  distance  is 
normalized  by  the  radius  of  the  central  jet  (R)  and  the  axial  distance  is  normalized  by  the  central 
jet  diameter  (D). 

Measurements  [7-4]  were  taken  across  the  diameter  of  the  flow.  However,  in  the  figures 
presented  here,  the  data  on  either  side  of  the  centerline  are  shown  on  the  same  side.  This  gives 
an  indication  of  the  symmetry  of  the  flow  and  the  amount  of  scatter  in  the  data.  Measurements 
have  been  reported  for  x/D  =  0.16  (initial  plane),  1.06, 2.65, 5.29, 7.94, 15.9,  and  26.5.  Figure  7- 
1  shows  the  conditional  mean  axial  velocities  from  the  stochastic  model  at  four  different  axial 
stations.  The  calculations  are  in  very  good  agreement  with  data.  The  conditional  velocities  show 
the  momentum  transfer  between  the  jets  as  they  mix.  The  annular  jet  axial  velocity  is  accelerated 
as  the  annular  fluid  moves  towards  the  centerline  and  the  central  jet  decelerates.  Similarly,  the 
coflow  fluid  accelerates  as  it  moves  in.  The  conditional  velocities  differ  significantly  from  each 
other  in  the  developing  region,  and  these  differences  become  insignificant  far  downstream  (x/D  = 
26.5). 

The  corresponding  mean  axial  velocity  profiles  from  the  <to>-model  are  shown  in  Figure  7-2. 
The  profiles  from  the  <co>-model  show  greater  decay  of  the  centerline  velocity  than  those  shown 
by  the  stochastic  model  and  data.  However,  the  conditional  velocities  show  the  same  relative 
trends  and  the  agreement  with  data  over  most  of  the  radius  is  quite  good.  The  spreading  rate 
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Lines  -  JPDF  calculations 
Symbols  -  Experimental  data 

D  Central  jet  fluid 
•  Swirler  fluid 
v  Coflow  fluid 


x/D  =  1 .06 


5.29 


<U>/<U>oc 


Figure  7-1.  Radial  profiles  of  conditional  mean  axial  velocity  from  pdf  (stochastic 
velocity-frequency  model)  calculations  compared  against  data  at  various  axial  stations. 
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x/D  =  1 .06  5.29 


<U>/<U>oc 


Figure  7-2.  Radial  profiles  of  conditional  mean  axial  velocity  from  pdf  (<co>-model)  calculations 
compared  against  data  at  various  axial  stations  (see  Figure  7-1  for  legend). 
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(and  the  centerline  velocity  decay)  can  be  adjusted  by  changing  the  values  of  Cj  and  C2. 
However,  tests  indicated  that  changing  the  constants  to  get  a  better  agreement  on  the  spreading 
rate  would  also  lower  the  turbulent  kinetic  energy  (i.e.,  <u2>,  <v2>,  and  <w2>)  and  results  in  a 
poorer  agreement  with  data  for  those  and  other  turbulent  correlations. 

The  profiles  of  <co>  obtained  from  the  stochastic  and  <0)>-model  are  shown  in  Figure  7-3.  As 
explained  earlier,  the  <(0>  profiles,  especially  the  peak  values,  from  the  stochastic  model  are 
considerably  lower  than  those  from  the  <co>-model  in  the  developing  region  (x/D  <  7.94).  In 
this  region  the  production  of  <co>  plays  a  major  role  due  to  the  large  mean  velocity  gradients. 
The  production  is  lower  in  the  stochastic  model  due  to  the  choice  of  the  constants.  Further 
downstream,  production  is  less  and  the  dissipation  of  <co>  plays  a  significant  role.  Since  the 
dissipation  is  lower  in  the  stochastic  model  than  in  the  <co>-model  it  can  be  seen  that  the  <(0> 
values  are  higher  for  the  stochastic  model  at  x/D  =  26.5.  In  fact,  the  two  profiles  crossover  at 
about  x/D  =  16. 


x/D  =  1.06  5.29  7.94  26.5 


Figure  7-3.  Profiles  of  mean  relaxation  rate  <(o>  from  stochastic  model  calculations  (solid  lines) 
compared  with  those  from  <co>-model  calculations  (dashed  lines)  at  various  axial  stations. 
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The  profiles  of  the  mean  swirl  velocity  are  shown  in  Figure  7-4.  The  decay  of  the  swirl  velocity 
is  well  predicted  at  all  stations  including  the  downstream  stations  not  shown  here.  The  increase 
in  swirl  velocity  of  the  central  jet  and  coflow  fluids  as  they  move  towards  the  annular  swirling  jet 
is  clearly  depicted  by  both  the  measured  and  calculated  profiles.  The  profiles  of  <W>  calculated 
by  the  <co>-model  are  also  in  good  agreement  with  data.  The  negative  values  of  <W>  shown  by 
data  at  x/D  =  5.29  and  7.94  are  nonphysical  and  indicative  of  some  experimental  error 

Since  most  significant  differences  among  the  conditional  quantities  exist  in  the  near  nozzle 
region,  x/D  <  7.94,  only  selected  stations  in  this  region  are  shown,  in  the  interest  of  clarity,  in  the 
figures  to  follow.  The  results  are  in  good  agreement  with  data  at  all  stations.  The  profiles  of 
turbulent  kinetic  energy  from  the  stochastic  model  are  shown  in  Figure  7-5  for  x/D  =  5.29  and 
7.94. 


x/D  =  1 .06 


5.29 


7.94 


1 0<W>/<U>oc 


Figure  7-4.  Profiles  of  conditional  mean  swirl  velocity  <W>  from  stochastic  model  calculations 

compared  against  data  (see  Figure  7-1  for  legend). 


7-10 


x/D  =  5.29  7.94 


1 00k/<U>2oc 

Figure  7-5.  Profiles  of  conditional  turbulent  kinetic  energy  from  stochastic  model  calculations 
compared  against  data  (see  Figure  7-1  for  legend). 


Figure  7-6  shows  the  unconditional  and  conditional  mean  axial  velocity  and  turbulent  kinetic 
energy  at  x/D  =  5.29.  The  unconditional  mean  velocity  lies  within  the  conditional  velocities  at 
any  given  radius  since: 

<U>  =  Yj  <U>j  +  ya  <U>a  +  Yc  <U>C  (7.5) 

where  the  subscripts  j,  a,  and  c  denote  the  central  jet,  annular  jet,  and  coflow,  respectively,  and  y 
is  the  fraction  of  the  time  that  the  fluid  from  the  particular  jet  is  found  at  the  measurement 
location  (note  that  y's  are  less  than  unity  and  add  up  to  unity).  However,  the  unconditional 
turbulent  correlations,  in  general,  need  not  lie  within  the  conditional  quantities.  Expressions 
similar  to  Eq.  7.5  can  be  derived  for  each  of  the  correlations  that  will  have,  in  addition  to  a  linear 
combination  of  the  conditional  correlations,  terms  involving  differences  between  the  conditional 
mean  velocities  and  other  conditional  quantities.  It  is  easier  to  visualize  that  while  the 
conditional  fluctuations  are  calculated  from  the  corresponding  conditional  mean,  the 
unconditional  fluctuations  are  with  respect  to  the  unconditional  mean.  For  example,  in  a  region 
where  the  conditional  velocities  are  significantly  different,  the  root  mean  square  (rms)  fluctuation 
about  the  unconditional  mean  can  be  significantly  larger  than  the  rms  fluctuations  about  the 
respective  conditional  means.  This  fact  is  illustrated  by  the  kinetic  energy  profiles  in  Figure  7-6. 
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Figure  7-6.  Profiles  of  conditional  and  unconditional  mean  axial  velocity  and  turbulent  kinetic 
energy  at  x/D  =  5.29  from  stochastic  model  calculations.  Solid  line  -  jet  fluid;  short  dashed  line 
-  annular  jet  fluid;  long  dashed  line  -  coflow;  extra  long  dashed  line  -  unconditional. 


The  kinetic  energy  profiles  from  the  <co>-model  (not  shown  here)  are  also  in  good  agreement 
with  data,  but  some  differences  in  the  details  exist  between  the  results  from  the  two  models  as 
illustrated  by  the  next  two  figures.  Figures  7-7  and  7-8  show  the  conditional  axial  velocity 
variance  from  the  stochastic  and  <to>-models  respectively.  Both  models  produce  results  that 
have  the  same  level  of  agreement  with  data.  However,  the  shift  in  the  location  of  the  peaks  for 
the  central  and  annular  jets  is  better  captured  by  the  stochastic  model  than  the  <co>-model.  This 
may  be  due  to  the  fact  that  the  turbulent  transport  of  CD  is  modeled  (by  gradient  diffusion)  in  the 
<co>-model  while  it  need  not  be  modeled  with  the  stochastic  model. 

This  discrepancy  in  the  <u2>  profile  is  further  amplified  in  the  <u^>  profiles  shown  in  Figures  9 
and  10.  The  quantity  <u2v>  can  be  interpreted  as  the  radial  turbulent  flux  of  <u2>.  The  quantity 
<u2v>  will  have  opposite  signs  on  either  side  of  the  location  (approximately)  of  the  peak  of  <u2> 
since  turbulence  will  disperse  the  peak  in  opposite  directions  from  the  peak.  Figure  7-9  shows 
that  the  stochastic  model  shows  good  agreement  with  data  while  the  <co>-model  calculations  are 
not  in  as  good  an  agreement  (Figure  7-10). 
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Figure  7-7.  Profiles  of  conditional  axial  velocity  variance  from  stochastic  model  calculations 
compared  against  data  (see  Figure  7-1  for  legend).. 
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Figure  7-8.  Profiles  of  conditional  axial  velocity  variance  from  <to>-model  calculations 
compared  against  data  (see  Figure  7-1  for  legend).. 
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Figure  7-9.  Profiles  of  (conditional)  triple  correlation  <u2v>  from  stochastic  model  calculations 

compared  against  data  (see  Figure  7-1  for  legend). 
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Figure  7-10.  Profiles  of  (conditional)  triple  correlation  <u2v>  from  <co>-model  calculations 
compared  against  data  (see  Figure  7-1  for  legend). 
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One  can  draw  an  analogy  between  <u  v>  (the  transport  of  the  velocity  variance  <u  >  in  the 

v-direction)  and  «|)2v>  (the  transport  of  <(|)2>,  variance  of  a  scalar  <J)  in  the  v-direction).  It  is 

important  to  calculate  scalar  transport  (both  mean  and  fluctuation  of  the  scalar)  accurately, 

especially  in  reacting  flows,  since  it  significantly  affects  local  heat  release  rates  and  pollutant 

formation.  It  is  even  more  critical  when  there  is  mixing  between  different  streams  (such  as  fuel 

and  oxidizer  in  diffusion  flames). 

9 

In  Reynolds-stress  closures,  all  third-order  correlations  (such  as  <u  v>  )  are  modeled,  usually  by 
gradient  diffusion.  In  the  present  study,  the  flow  under  consideration  was  calculated  using  the 
RSM  of  Daly  and  Harlow  [7-7]  that  is  widely  used.  The  details  of  the  code  and  models  used  are 
presented  in  Ref.  7-8.  The  model  for  <u  v>  is: 


<u2v> 


-C. 


k 

<£> 


<v2> 


9<U> 

9r 


(7.6) 


where  Cs  =  0.22.  The  profile  of  (unconditional)  <u2v>  calculated  from  Eq.  7.6  with  the 
right-hand  side  quantities  taken  from  the  Reynolds-stress  calculations  (short  dashed  line)  is 
compared  against  the  unconditional  <u2v>  profile  from  the  stochastic  model  calculations  (solid 
line)  in  Figure  7-11.  The  comparison  shows  that  the  RSM  calculations  severely  underpredict  the 
transport  of  <u2>  especially  at  the  stations  close  to  the  nozzle.  The  magnitude  of  <u2v>  is  lower 
because  the  RSM  calculations  also  grossly  underpredict  <u  >  in  this  region.  In  addition  the 
<u2v>  profiles  from  the  RSM  calculations  change  sign  at  a  different  radial  location  than  the  pdf 
calculations  since  the  peaks  of  <u2>  are  also  shifted  in  the  RSM  calculations.  These  differences 
get  larger  further  downstream.  A  better  test  of  the  gradient  diffusion  hypothesis  would  be  to  use 
the  pdf  results  for  the  right-hand  side  in  Eq.  7.6.  These  profiles  (long  dashed  lines  in  Figure  7- 
11)  also  show  differences  both  in  magnitude  and  the  location  where  the  profiles  change  sign. 
Though  these  discrepancies  are  not  dramatic  for  the  flow  considered  here,  they  underscore  the 
fact  that  the  gradient  diffusion  hypothesis  could  lead  to  different  results  and  in  the  case  of 
reacting  flows  where  counter-gradient  diffusion  has  been  demonstrated  the  results  would  be 
erroneous. 


It  should  be  noted  that  even  with  the  scalar  pdf  method  in  which  the  joint  pdf  of  scalars  only  is 
considered,  the  turbulent  flux  of  the  scalar  pdf  is  modeled  by  gradient-diffusion.  While  the 
scalar-pdf  method  is  a  useful  tool  for  the  development  of  reduced  reaction  mechanisms,  since 
reaction  still  appears  in  closed  form  in  the  method,  the  full  potential  of  the  pdf  method  can  only 
be  realized  by  the  velocity- scalar  joint  pdf  method. 
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Figure  7-11.  Profiles  of  (unconditional)  triple  correlation  <u2v>  from  stochastic  model  calculations 
(solid  line);  from  Reynolds-stress  model  (Eq.  7.6)  using  results  from  Reynolds-stress  closure 
calculations  (short  dashed  line);  and  from  model  (Eq.  7.6)  using  results  from  stochastic  model 

calculations  (long  dashed  line). 


Figure  7-12  shows  the  conditional  <u3>  profiles  at  x/D  =  2.65  and  5.29.  The  dramatic 
differences  between  the  conditional  quantities  is  well  predicted  by  the  calculations.  The  results 
from  the  <co>-model  show  similar  discrepancies  as  seen  for  <u  v>.  Figure  7-13  shows  the 
(conditional)  fourth  order  correlation  <w4>  from  the  stochastic  model  at  x/D  =  5.29  and  7.94. 
Considering  that  fourth  or-der  conditional  statistics  are  being  measured  and  computed,  the  results 
are  in  very  good  agreement  with  data.  The  corresponding  profiles  for  the  <co>-model  are  shown 
in  Figure  7-14.  These  profiles  are  also  in  good  agreement  with  data.  However  the  radial 
separation  between  the  peak  locations  for  the  central  jet  and  annulus  at  x/D  =  5.29  are  better 
predicted  by  the  stochastic  model. 
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Figure  7-12.  Profiles  of  conditional  <u3>  from  stochastic  model  compared  against  data  (see 

Figure  7-1  for  legend) 
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Figure  7-13.  Profiles  of  (conditional)  fourth-order  correlation  <w4>  from  stochastic  model 
calculations  compared  against  data  (see  Figure  7-1  for  legend). 
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Figure  7- 14.  Profiles  of  (conditional)  fourth  order  correlation  <w4>  from  <co>-model  compared 

against  data  (see  Figure  7-1  for  legend). 


7.6  Summarizing  Remarks 

Swirling  flow  calculations  have  been  made  with  the  pdf  method.  The  newly  developed 
stochastic  frequency  and  velocity  models  (SFM  I  and  RLM  respectively)  have  been  validated 
against  benchmark  quality  data  for  a  coaxial  swirling  jet  flow.  Detailed  comparisons  between 
measured  and  calculated  conditional  means  and  correlations  up  to  fourth  order  have  been  made. 
The  calculations  are  in  very  good  agreement  with  data.  The  present  study  suggests  the  values 
C(oi  =  0.04  and  Cffi2  =  0.7  for  the  two  constants  that  significantly  affect  the  evolution  of  the  mean 
frequency  <co>  in  the  SFM  I  model. 

Swirling  flows  were  calculated  for  the  first  time  using  the  pdf  method.  While  extensions  to  the 
algorithm  were  made  to  include  the  axial  pressure  gradient  and  additional  terms  due  to  nonzero 
mean  swirl  velocity,  no  specific  modifications  were  made  to  the  SFM  and  RLM  models  or 
elsewhere  on  account  of  swirl. 

Calculations  were  also  made  with  the  pdf  method  using  a  transport  model  for  <co>.  These  results 
are  also  in  reasonably  good  agreement  with  data,  although  the  stochastic  model  better  predicts 
some  of  the  details  of  turbulent  transport. 
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The  calculated  results  from  the  stochastic  model  were  compared  against  gradient-diffusion 
models  used  by  conventional  Reynolds-stress  closures,  and  the  latter  were  found  to  be  deficient. 
By  extension,  the  scalar  pdf  method  is  also  subject  to  the  same  limitations  since  gradient 
diffusion  modeling  is  used  for  turbulent  transport  of  the  scalar  pdf. 

The  stochastic  frequency  model  provides  complete  closure  for  the  joint  pdf  method  and 
eliminates  the  need  to  supply  a  time  scale  for  the  solution  of  the  pdf  transport  equation.  The 
stochastic  model  also  contains  more  information  than  that  based  on  <co>,  and  allows  the 
flexibility  to  build  in  more  information  than  currently  used  in  further  developing  the  model.  The 
application  of  the  SFM  I  and  RLM  models  to  swirling  jet  diffusion  flames  is  reported  in  the  next 
section. 
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8.  SWIRLING  JET  DIFFUSION  FLAMES 


8.1  Overview 


Computations  using  the  joint  velocity-frequency-scalar  pdf  method  as  well  as  benchmark  quality 
experimental  data  for  swirling  and  nonswirling  hydrogen  jet  diffusion  flames  are  reported  in  this 
chapter.  Previous  studies  of  diffusion  flames  reported  in  literature  have  been  limited  to 
nonswirling  flames  and  have  had  no  detailed  velocity  data  reported  in  the  developing  (near¬ 
nozzle)  region  of  the  flames.  The  measurements  and  computations  reported  herein  include 
velocities  (mean  and  higher  moments  up  to  fourth  order)  and  temperature  (mean  and  variance) 
near  the  burner  exit  and  downstream  locations  up  to  26.5  jet  diameters.  The  computed  results  are 
in  good  agreement  with  data.  This  study  serves  to  further  validate  the  joint  velocity-ffequency- 
composition  pdf  method  for  reacting  flows  as  well  as  present  data  that  can  be  used  for  validation 
of  other  reacting  flow  models  as  well. 


8.2  Background 

Although  jet  diffusion  flames  have  been  studied  extensively  and  reported  in  the  literature,  there 
are  voids  with  respect  to  two  key  aspects.  First,  there  are  no  detailed  measurements,  especially 
velocity  measurements,  in  the  developing  regions  of  the  flames,  and  second,  the  flames  have 
been  predominantly  nonswirling.  Swirl  is  widely  used  in  practical  combustion  systems  such  as 
gas  turbine  combustors  for  flame  stabilization  and  enhancing  fuel-air  mixing  and  combustion 
intensity.  It  is  necessary  that  the  turbulent  combustion  models  be  validated  in  detail  in  simpler 
flames  incorporating  the  essential  features  of  practical  flows  before  they  can  be  used  with 
confidence  to  design  practical  combustion  systems.  In  the  present  study  detailed  velocity  and 
temperature  measurements  are  presented  for  swirling  and  nonswirling  hydrogen  jet  diffusion 
flames  in  the  developing  region  of  the  flames  (<30  jet  diameters  in  axial  distance),  and  the 
measurements  are  used  to  further  validate  the  joint  pdf  method. 

The  present  computations  of  diffusions  flames  with  the  pdf  method  differ  from  previous 
computations  in  that  for  the  first  time  a)  swirling  diffusion  flames  are  being  computed  and  b) 
detailed  comparison  with  data  is  made  in  the  developing  region  of  the  flame.  As  explaied  in  the 
previous  chapter,  notable  feature  of  the  measurements  reported  in  this  study  is  that  the  velocity 
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correlations  are  conditional  upon  the  fluid  originating  from  a  given  inlet  stream.  As  a  result 
errors  due  to  velocity  bias  are  avoided.  These  conditional  quantities  are  readily  computed  in  the 
pdf  method  without  need  for  additional  modeling,  and  serve  as  severe  tests  for  the  method’s 
ability  to  predict  the  details  of  the  individual  streams  and  the  transport  processes  between  the 
streams. 

The  experiments  were  conducted  by  Takahashi  et  al.  (UDRI)  at  WPAFB’s  Fuels  Laboratory. 
The  study  is  reported  in  detail  in  Ref.  8-1. 


8.3  Experimental  Setup  and  Techniques 

Figure  8-1  shows  a  schematic  of  the  combustor  system  used.  It  consists  of  a  central  fuel  tube 
(9.45  mm  inner  diameter  [D],  0.2  mm  lip  thickness,  806  mm  length)  and  a  concentric  annulus-air 
tube  (26.92  mm  inner  diameter)  centered  in  a  vertical  test  section  (150  x  150  mm  square  cross 
section  with  rounded  corners,  486  mm  length)  through  which  external  air  is  supplied.  The  test 
section  is  sided  with  four  quartz  windows  for  optical  observations  and  diagnostics.  A  helical 
vane  swirler  unit  can  be  placed  in  the  annulus  channel  96  mm  upstream  from  the  jet  exit. 
Swirlers  of  vane  angles  0  deg,  30  deg,  and  45  deg  were  used  in  the  present  study. 

The  three-component  LDV  system,  described  in  detail  in  Takahashi  et  al.  [8-2],  consists  of  the 
two  segments  of  three-beam  two-channel  optics  and  two-beam  one-channel  optics.  The  former 
utilizes  a  514.5  nm  line  of  an  argon-ion  laser  (Spectra  Physics  171;  15W  nominal  output)  and 
measures  the  velocity  components  in  the  directions  of  ±  45  deg  off  the  jet  axis.  The  latter 
utilizes  a  488.0  nm  line  of  the  laser  and  measures  the  tangential  velocity  component.  The 
overlapping  probe  volume  is  approximately  a  100  (xm  diameter  sphere.  The  calculated  fringe 
spacing  is  approximately  3.6  mm.  Submicron-size  zirconia  (ZrO^  particles  (97%,  <1  jxm)  are 
used  as  the  seed  particles.  A  portion  of  the  data  was  filtered  out  by  the  so-called  n-G  method 
(i.e.,  velocities  whose  deviation  from  the  mean  exceeded  n  times  the  standard  deviation  [g]  are 
eliminated).  The  coefficient  n  =  4  was  employed  because  preliminary  processing  tests  showed 
that  n  =  3  cut  off  some  valid  data  and  altered  high  moments  considerably. 

LDV  measurements  were  made  by  seeding  one  stream  at  a  time  and  completing  all  the  axial  and 
radical  scans  before  seeding  the  next  stream  and  repeating  the  scans.  All  the  LDV  measurements 
reported  are  conditional  upon  the  fluid  originating  from  either  the  jet  annulus  or  the  coflow. 

The  CARS  system,  described  in  detail  in  Takahashi  et  al.  [8-3],  consists  of  a  pulsed  NdrYAG 
laser  (Quanta  Ray  DCR-2A,  10  ns  pulse  width,  10  FIz  repetition  rate),  dye  laser  optics,  incident 
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Figure  8-1.  Schematic  of  the  Swirling  Jet  Diffusion  Flame  Combustor  test  section. 


and  collection  optics,  a  3/4  m  grating  spectrometer,  and  an  intensified  charge-coupled  device 
(CCD)  camera  (Princeton  Instruments).  The  output  from  the  laser  is  frequency-doubled  (532 
nm,  ~150  mJ)  and  divided  into  four  beams  of  nearly  equal  intensity.  Two  of  these  serve  as  the 
pump  beams,  while  the  other  two  pump  a  dye  laser  oscillator  and  amplifier  to  provide  a 
broadband  Stokes  beam  centered  at  607  nm.  The  Stokes  beam  and  the  two  pump  beams  are  then 
focused  together  in  a  folded  BOXCARS  configuration.  The  effective  probe  volume  size  is 
estimated  at  approximately  25  mm  in  diameter  and  250  mm  in  length.  Typically,  500  CARS 
signals  are  acquired  at  each  location  and  processed  by  a  microcomputer.  The  CARS 
measurements  are  based  on  nitrogen  molecules  originated  in  air.  Unlike  seed  particles  which 
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follow  fluid  elements,  molecules  diffuse  among  different  species.  Therefore,  the  CARS  data  are 
unconditional  upon  the  origin  of  the  fluid  if  they  are  mixed.  In  the  jet-to-annulus  fluid  boundary 
zone  near  the  jet  exit,  the  determination  of  CARS  temperature  becomes  difficult  due  to  the 
interference  by  nonresonant  background  emission  from  hydrogen  molecules.  The  accuracy  of 
the  temperature  measurements  is  estimated  to  range  between  10%  near  room  temperature  and 
5%  near  the  flame  temperature,  with  the  largest  contribution  to  uncertainty  from  shot-to-shot 
variation  in  the  Stokes-laser  spectral  distribution. 


8.4  Computations 

8.4.1  The  Joint  PDF  Method  and  Models  Used 

The  pdf  method  used  in  the  present  study  is  the  same  as  that  used  for  the  constant-density 
swirling  flow  study  reported  inthe  previous  chapter,  except  for  the  addition  of  the 
thermochemistry  and  molecular  diffusion  due  to  mean  scalar  gradient  described  later  and, 
obviously,  the  variability  of  density.  The  SFM  I  and  RLM  models  are  used  for  the  frequency 
and  velocity  respectively  (see  sections  3.3  and  3.4). 

The  boundary-layer  (parabolic)  algorithm  is  used  here.  At  the  axial  position  x,  each 
computational  particle  has  the  position  x*(x),  velocity  U*(x),  the  relaxation  rate  or  frequency 
oo*(x),  and  four  scalar  values  $*(x),  p*(<t>*),  T*(<|>*)  and  c*(x).  The  scalar  <J>*  is  the  mixture 
fraction  (conserved  scalar)  used  to  model  the  thermochemistry.  Its  value  changes  due  to 
turbulent  and  molecular  mixing  and  molecular  diffusion.  The  scalars  p*  and  T*  are  density  and 
temperature,  respectively,  and  are  derived  from  0*  based  on  the  thermochemistry  used.  The 
scalar  c*  is  a  passive  scalar  with  no  source  terms  (i.e.,  c*  is  held  constant),  and  is  used  to  tag 
particles  by  their  jet  of  origin. 

The  turbulent  mixing  (i.e.  turbulent  convection)  of  the  mixture  fraction  (<(>*)  is  in  closed  form  and 
need  not  be  modeled.  The  molecular  mixing  is  modeled  using  the  improved  mixing  model 
(IMM)  with  the  standard  model  constant  =  2.0.  The  molecular  diffusion  process  is 
represented  by  Fick  s  Law.  Typically,  although  it  is  in  closed  form,  molecular  diffusion  due  to 
mean  scalar  gradient  is  not  explicitly  included  in  the  computations  for  high  Reynolds  number 
turbulent  flows  since  turbulent  mixing  is  dominant.  However,  in  the  present  study,  molecular 
diffusion  is  explicitly  included  since  the  stoichiometric  surface  lies  at  the  outer  edge  of  the  shear 
layer  (for  most  diffusion  flames  and  especially  for  hydrogen  flames  due  to  the  low  value  of  the 
stoichiometric  mixture  fraction  and  the  high  diffusivity  of  hydrogen).  In  this  region,  the 
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turbulent  intensity  is  low  and  molecular  diffusion  plays  an  important  role  in  the  radial  spread  of 
the  flame.  The  change  in  0*  over  the  particle  time-step  At*due  to  molecular  diffusion  is  given 
by: 


Atf)  * 


At*  d 
p*  8r 


(' 


<p>D, 


9<0>  \ 
dr  )  , 


(8.1) 


where  <p>  and  <0>  are  the  mean  density  and  mixture  fraction  respectively,  and  is  the 
diffusivity  taken  as  a  function  of  <0>  from  Bilger  [8-4]: 

D<j>  =  0.26  +  257  «j»  +  3.5  x  108  <(>>5  (cm2/s)  for  0  <  «J»  <  0.0325;  and 

D<j,  =  3.8  -  3.03  «]»  +  27.5  (1.0001  -  <0>)10  (cm2/s)  for  0.0325  <  <0>  <  1.0. 

Means  and  other  correlations  are  extracted  from  the  solution  at  any  desired  axial  station  x  by 
forming  sums  of  particle  properties  within  spatial  bins  in  the  lateral  direction  (r-direction)  and 
fitting  cubic  B-splines  to  these  sums.  Typically  density-weighted  (Favre)  means  are  calculated 
as  is  appropriate  for  reacting  flows;  however,  means  weighted  by  any  desired  quantity  can  be 
computed. 


8.4.2  Thermochemistry 

Due  to  the  highly  reactive  nature  of  the  fuel  used,  namely  hydrogen,  the  thermochemistry  in  the 
present  study  is  described  using  a  fast  one- step  reaction  in  equilibrium.  This  allows  the 
modeling  of  the  diffusion  flame  using  a  single  conserved  scalar,  the  mixture  fraction  which  is 
unity  in  the  hydrogen  stream  and  zero  in  the  air  streams.  The  reaction  is  completely  described  by 
the  evolution  of  the  mixture  fraction.  The  temperature  and  density  are  then  functions  of  mixture 
fraction  alone  and  are  determined  from  the  equilibrium  composition  at  the  given  value  of  <|>, 
using  the  heat  of  combustion,  the  specific  heat  at  constant  pressure  (Cp)  of  the  mixture  as  a 
function  of  temperature,  and  the  ideal  gas  law.  In  the  present  study  the  equilibrium  calculations 
were  performed  at  different  values  of  (|)  prior  to  the  pdf  calculations  and  the  values  of  specific 
volume  (inverse  of  density)  and  temperature  were  tabulated  for  use  (through  linear  interpolation) 
in  the  pdf  calculations.  The  tables  consisted  of  21  points  between  0  =  0  and  the  stoichiometric 
value  0st  =  0.0282  and  51  points  between  0st  and  0  =  1  to  provide  sufficient  accuracy  for  the  pdf 
calculations. 
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8.4.3  Initial  Conditions 


Initial  conditions  for  the  computations  are  prescribed  from  experimental  data.  The  first 
measurement  station  x  =  1.5  mm  (x/D  =0.16)  is  taken  as  the  initial  plane  for  the  calculations. 
The  initial  velocity  pdf  is  prescribed  to  be  joint  normal  with  the  mean  and  covariances  taken 
from  linearly  interpolated  experimental  data.  The  initial  pdf  of  turbulence  frequency  is  taken  to 
be  log-normal,  i.e.,  ln(co*/<co>)  has  a  normal  distribution  with  mean  and  variance  -1/2  and  1, 
respectively,  in  accordance  with  the  construction  of  the  model.  The  initial  profile  of  <co>  is 
derived  from  experimental  data  using  the  following  expression,  as  described  in  the  previous 
chapter. 

The  initial  profiles  of  the  mean  mixture  fraction  and  its  variance  «|),2>  are  deduced  from  the 
measured  mean  and  variance  of  temperature  at  the  x  =  1.5  mm  location.  The  value  of  «j)>  is 
deduced  from  the  tabulated  values  of  temperature  versus  mixture  fraction  and  the  value  of  <(j)’  > 
is  estimated  by 


<§'2> 


(8.2) 


where  T'2  is  the  measured  variance  of  temperature,  and  Tst  (=  2377K)  is  the  temperature  at  the 
stoichiometric  mixture  fraction  (J)st.  The  pdf  of  <j)  is  prescribed  to  be  Gaussian  with  <())>  and 
variance  «|>  >. 


8.5  Results  and  Discussion 


Experiments  and  computations  were  performed  for  the  four  cases  listed  in  Table  8-1.  The 
velocities  Uj,  Ua,  and  Ue  are  the  nominal  bulk-averaged  axial  velocities  of  the  fuel  jet,  annular 
air  jet,  and  the  coflowing  air  stream,  respectively.  Two  of  the  cases  are  nonswirling  and  two  are 
swirling.  For  Case  2,  the  high-velocity  nonswirling  case,  the  lip  thickness  of  the  fuel 
tube  used  was  1.2  mm  in  order  to  anchor  the  flame  and  keep  it  from  lifting.  For  all  other  cases, 
the  tube  walls  were  chamfered  to  form  knife  edges  at  the  exit 
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TABLE  8-1. 
FLOW  CONDITIONS 


Case  No. 

Uj  -  m/s 

Ua  -  m/s 

Uo  -  m/s 

Swirler  Angle 

1 

25 

4 

1 

0 

2 

100 

20 

4 

0 

3 

100 

20 

4 

30 

4 

100 

20 

4 

45 

The  nominal  Reynolds  number  for  the  hydrogen  jet  based  on  Uj  and  the  jet  diameter  (D)  is  2208 
for  Case  1  and  8830  for  the  others,  and  the  nominal  Reynolds  number  for  the  annular  air  jet 
based  on  Ua  and  its  hydraulic  diameter  is  8400  for  Case  1  and  33,600  for  the  others.  The  swirl 
numbers  for  the  annular  jet,  calculated  from  the  measured  axial  and  tangential  velocities  at  the 
initial  plane,  are  0.382  for  Case  3  and  0.516  for  Case  4. 

To  be  concise,  only  sample  results  from  a  nonswirling  case  (Case  2)  and  the  two  swirling  cases 
in  Table  8-1  will  be  presented  here.  Case  2  was  selected  to  allow  a  more  direct  comparison 
between  the  nonswirling  and  swirling  cases.  Detailed  tables  and  plots  of  the  measured  data  for 
the  four  cases  are  reported  in  Takahashi  et  al.  [8-2.  8-3,  8-5  through  8-8],  and  the  tabulated  data 
are  available  on  diskettes.  Measurements  have  been  made  at  axial  stations  x/D  =  0.16  (initial 
plane),  1.06,  2.65,  5.29,  7.94, 15.9,  and  26.5.  Radial  profiles  of  conditional  mean  velocities  and 
velocity  correlations  up  to  fourth  order  as  well  as  unconditional  mean  and  variance  of 
temperature  have  been  reported. 

The  notation  used  to  present  the  results  for  the  velocity  components  is  U  =  (U,  V,  W)  =  (<U>  + 
u,  <V>  +  v,  <W>  +  w).  The  notation  for  the  temperature  results  are  presented  later.  The 
measured  mean  axial  velocity  on  the  centerline  at  x/D  =  0.16  denoted  as  <U>0C  is  used  to 
normalize  the  velocities  and  the  temperature  at  the  stoichiometric  mixture  fraction,  Tst  (= 
2377K),  is  used  to  normalize  temperatures.  The  values  of  <U>0C  are  132.8  m/s,  130.3  m/s,  and 
130.0  m/s  for  Cases  2,  3,  and  4,  respectively.  The  radial  distance  is  normalized  by  the  fuel  jet 
radius  R  (=  D/2). 

Figure  8-2  shows  the  computed  conditional  and  unconditional  mean  axial  velocity  profiles  for 
the  nonswirling  case  (Case  2)  compared  against  measurements  at  various  axial  stations.  The 
calculations  are  in  good  agreement  with  data.  There  are  differences  in  the  velocities  of  the  fluids 
originating  from  the  different  streams  as  expected.  The  annulus  air  increases  its  axial  velocity  as 
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it  moves  towards  the  central  jet,  but  its  axial  velocity  remains  lower  than  that  of  the  central  jet. 
As  the  jets  spread,  the  profiles  become  flatter  and  the  differences  get  smaller  at  the  downstream 
locations  x/D  =  15.9  and  x/D  =  26.5  (not  shown).  While  the  unconditional  data  cannot  be 
deduced  from  the  conditional  data  without  additional  measurement  of  the  intermittency  factor  for 
the  different  streams,  unconditional  and  conditional  quantities  of  any  order  can  be  extracted  from 
the  pdf  solution.  The  unconditional  mean  axial  velocity  in  Figure  8-2  lies  between  the 
conditional  values  as  expected. 

Figure  8-3  shows  the  computed  profiles  of  mean  temperature  for  Case  2  compared  against  data. 
Both  the  Favre-averaged  (density-weighted)  mean  temperature  <T>  as  well  as  the  Reynolds- 

averaged  (volume-weighted  or  spatially  averaged)  mean  temperature  denoted  by  T  are  shown. 
The  Favre  average  is  lower  in  the  flames  outer  radial  region,  where  there  is  significant 
probability  of  hot  products  and  cold  air  since  the  cold  air  is  an  order  of  magnitude  denser  (1.17 
kg/m3  at  the  inlet  temperature  of  298K)  and  contributes  more  to  the  average  than  the  hot  rarer  air 
(0.125  kg/m3  at  Tst),  whereas  on  the  inside  of  the  flame  the  difference  is  negligible  since  the 
densities  of  cold  unburnt  hydrogen  (0.082  kg/m3  at  the  inlet  temperature  of  298K)  and  the  burnt 
products  are  of  the  same  order  of  magnitude.  Clearly,  the  spatial  mean  is  much  closer  to  the 
measured  data,  illustrating  that  the  CARS  measurement  represents  a  spatial  average  and  is  not 
weighted  by  the  mass  flow  into  the  measurement  volume.  This  is  an  important  observation  and 
draws  attention  to  the  fact  that  although  Favre  averages  are  conventionally  used  in  the 
computations  of  variable-density  flows,  one  has  to  be  careful  to  use  the  appropriate  averaging 
while  comparing  with  measurements.  On  the  other  hand,  if  the  seeding  of  the  flow  is  uniform, 
LDV  measurements  represent  Favre  means,  and  hence  the  Favre  means  are  used  while 
comparing  velocity  data.  However,  the  velocity  statistics  computed  using  Reynolds  averaging 
were  nearly  identical  to  those  computed  using  Favre-averaging.  This  indicates  that  the  velocity 
distribution  of  the  high  and  low  density  fluids  in  the  mixture  are  nearly  identical  at  any  point  in 
the  flow. 


The  calculated  profiles  of  (volume  averaged)  temperature  in  Figure  8-3  show  the  locations  of  the 
temperature  peaks  are  well  predicted;  however,  the  peak  values  are  lower  than  the  measured 
values  in  the  region  between  x/D  =  2.65  and  x/D  =  7.94.  At  the  further  downstream  locations 
x/D  =  15.9  and  x/D  =  26.5,  the  calculated  temperature  is  higher  than  the  data  at  the  outer  radial 
locations. 

Figure  8-4  shows  the  temperature  variance  for  the  nonswirling  case  (Case  2).  Again,  the 
Reynolds-averaged  variance  is  in  much  better  agreement  with  the  data  than  the  Favre  average. 
The  second  peak  of  temperature  variance  on  the  inside  (i.e.,  towards  the  centerline)  seen  in  the 
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calculations  corresponds  to  the  inner  gradient  of  the  temperature  peak.  Obviously,  this  is  the 
region  which  lacks  nitrogen  and  the  CARS  data  are  either  not  available  or  the  measured 
temperature  variance  is  expected  to  be  low  in  this  region. 

The  data  as  well  as  the  computations  in  Figures  8-2,  8-3,  and  8-4  show  that  the  flame  lies  at  the 
outside  edge  of  the  shear  layer  between  the  fuel  jet  and  the  annulus  jet  as  expected  due  to  the  low 
value  of  <|)st. 

Results  for  the  30  deg  swirl  case  (Case  3)  will  now  be  presented.  Figure  8-5  shows  the  radial 
profiles  of  mean  conditional  axial  velocity  at  different  downstream  stations  in  the  developing 
region.  For  clarity  of  the  figure,  the  calculated  unconditional  velocity  is  omitted.  The  spreading 
of  the  jet  is  much  faster  than  in  Case  2  as  a  result  of  the  enhanced  mixing  due  to  swirl  in  Case  3. 
The  location  of  the  shear  layer  is  at  a  larger  radial  distance  and  the  width  of  the  shear  layer  is 
also  larger  than  for  the  nonswirling  case  (Figure  8-2),  for  example  at  x/D  =  7.94.  Differences 
between  the  conditional  means  from  the  three  streams  are  evident. 

Overall,  the  computations  are  in  good  agreement  with  the  data  in  terms  of  the  trends,  locations, 
and  magnitudes  of  the  conditional  velocities;  however,  some  differences  exist.  The 
computations  show  a  slightly  slower  rate  of  spreading  of  the  hydrogen  jet  than  the  data. 
Although  the  spreading  rate  in  the  computations  could  be  adjusted  by  adjusting  the  constants  in 
the  Langevin  equation  and  the  turbulent  frequency  equation,  the  differences  are  not  significant 
enough  to  undertake  such  an  exercise;  hence,  the  constants  used  in  previous  study  (Chapter  7) 
were  retained.  The  data,  especially  at  x/D  =  2.65  and  5.29,  indicate  the  outer  edge  of  the 
hydrogen  jet  spreads  especially  faster  than  the  bulk  of  the  jet.  The  differences  between  the 
conditional  velocities  and  the  differences  between  the  computations  and  the  data  diminish  at  the 
downstream  locations  x/D  =  15.9  and  26.5  (not  shown  here). 

Figure  8-6  shows  the  conditional  mean  tangential  (or  swirl)  velocities  for  Case  2.  The 
development  and  decay  of  the  swirl  in  the  three  streams  is  well  predicted,  although  the 
computations  show  larger  differences  between  the  swirl  velocities  of  the  different  streams  than 
indicated  by  the  data.  The  computations  show  that  the  swirl  velocities  for  the  initially 
nonswirling  hydrogen  jet  and  coflowing  stream  lag  behind  those  for  the  annular  swirling  jet  as 
one  may  expect. 

The  calculated  Reynolds-averaged  mean  temperatures  T  are  compared  against  data  in  Figure  8-7. 
The  agreement  between  the  computations  and  the  data  is  better  than  in  Case  2  in  terms  of  the 
magnitude  and  shape  of  the  temperature  peaks.  For  the  swirling  case  (Case  3),  the  measured 
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Figure  8-2.  Radial  profiles  of  computed  Figure  8-3.  Computed  profiles  of 

conditional  and  unconditional  mean  axial  Favre-averaged  mean  temperature  (<T>)  and 
velocity  (lines)  compared  against  data  (symbols)  Reynolds-averaged  mean  temperature  (T) 
for  the  nonswirling  case  (Case  2).  compared  against  data  (symbols)  for  the 

nonswirling  case  (Case  2). 
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Figure  8-4.  Computed  profiles  of 
Favre-averaged  temperature  variance  (<T  >) 
and  Reynolds-averaged  temperature  variance 

(T'2)  compared  against  data  (symbols)  for  the 
nonswirling  case  (Case  2). 


Figure  8-5.  Radial  profiles  of  computed 
conditional  mean  axial  velocity  (lines) 
compared  against  data  (symbols)  for  the  30 
degree  swirl  case  (Case  3). 
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Figure  8-6.  Radial  profiles  of  computed 
conditional  mean  tangential  velocity  (lines) 
compared  against  data  (symbols)  for  the  30 
degree  swirl  case  (Case  3). 


Figure  8-7.  Computed  profiles  of 
Reynolds-averaged  mean  temperature  with 
(solid  lines)  and  without  (dotted  lines)  the 
inclusion  of  molecular  diffusion  compared 
against  data  (symbols)  for  the  30  degree  swirl 
case  (Case  3). 
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temperature  peaks  are  lower  at  corresponding  axial  stations  than  in  Case  2  (see  Figure  8-3)  and 
the  profiles  show  more  spreading  indicative  of  enhanced  turbulent  mixing  due  to  swirl.  Also,  at 
the  corresponding  axial  stations  the  peaks  are  located  at  a  slightly  larger  radial  distance  for  Case 
3  than  for  Case  2  due  to  the  centrifugal  action  of  the  swirl.  These  observations  are  consistent 
with  the  observations  made  for  the  mean  axial  velocity  profiles  shown  in  Figure  8-5. 

Figure  8-7  also  includes  the  calculated  profiles  of  T  without  the  explicit  inclusion  of  molecular 
diffusion  in  the  computations.  Overall,  the  inclusion  of  molecular  diffusion  improves  the 
agreement  of  the  computations  with  the  data.  It  is  interesting  to  note  that  the  computed  mean 
and  other  velocity  statistics  were  nearly  the  same  with  or  without  the  inclusion  of  molecular 
diffusion  although  the  mean  temperature  profiles  (and  other  scalar  statistics)  changed.  This 
indicates  that  small  changes  in  the  local  heat  release  rate  do  not  have  a  very  strong  influence  on 
the  hydrodynamics  of  the  flames  studied. 

The  measured  and  calculated  temperature  variances  are  shown  in  Figure  8-8.  The  comments 
regarding  the  inner  peak  in  the  computations  made  in  conjunction  with  Figure  8-4  are  more 
clearly  illustrated  in  this  figure.  The  magnitudes  and  shapes  of  the  outer  peaks  are  well 
predicted,  although  the  predictions  are  slightly  higher  for  the  downstream  locations  x/D  7.94. 

Figures  8-5  through  8-8  show  that  the  peak  temperature  (i.e.,  the  flame)  still  lies  at  the  outer  edge 
of  the  shear  layer  between  the  jet  and  the  annulus  flows,  but  the  broadening  of  the  temperature 
peaks  show  that  the  flame  is  spreading  into  the  shear  layer  more  than  in  the  nonswirling  case. 
Overall,  the  computations  and  data  (Figures  8-5  through  8-8)  show  enhanced  mixing  and 
spreading  of  the  flame  due  to  swirl  compared  to  Case  2.  The  computed  velocity  and  temperature 
statistics  are  in  good  agreement  with  the  data,  and  the  spreading  and  decay  of  the  mean  and  swirl 
velocities  are  well  predicted. 

Figures  8-9,  8-10,  and  8-11  show  the  conditional  mean  axial  and  swirl  velocities  and  the 
unconditional  Reynolds-averaged  temperature,  respectively  at  two  axial  stations  for  the  45  deg 
swirl  case  (Case  4).  The  axial  and  swirl  velocity  data  indicate  much  higher  mixing  and  spreading 
than  the  30  deg  swirl  case  (Case  3)  as  expected.  The  location  of  the  peak  temperature  (see  data 
in  Figure  8-11)  at  x/D  =  2.65  is  at  a  significantly  larger  distance  for  Case  4  (r/R  ~  2.5)  than  for 
Case  3  (r/R  ~  2.0),  while  at  x/D  =  7.94  the  peak  temperature  is  located  at  nearly  the  same 
location  (r/R  ~  2.15)  for  both  cases.  This  indicates  that  the  peak  flame  location  expands  to  a 
larger  radius  in  the  near  nozzle  region  for  the  45  deg  swirl  case,  but  moves  back  (indeed 
contracts)  to  about  the  same  location  as  for  the  30  deg  swirl  case,  at  the  downstream  locations. 
This  expansion  near  the  nozzle  is  due  to  the  larger  centrifugal  force  in  the  higher  swirl  case 
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forcing  both  the  annular  swirling  air  and  the  hydrogen  jet  to  a  larger  radius  (the  volume 
expansion  of  the  annulus  air  due  to  combustion  induces  a  radial  velocity  in  the  annular  jet  for  all 
the  cases).  Indeed,  the  data  for  mean  radial  velocities  at  the  initial  plane  (presented  in  Refs.  8-7 
and  8-8)  show  that  the  outer  edge  of  the  hydrogen  jet  and  the  bulk  of  the  annulus  air  have 
outward  mean  radial  velocities  of  approximately  5  m/s  for  Case  4  and  about  2  m/s  for  Case  3. 
Evidence  of  this  movement  can  be  seen  in  the  apparent  radial  shift  in  the  corresponding  radial 
profiles  of  the  mean  axial  and  tangential  velocities  at  x/D  =  2.65  in  Figures  8-9  and  8-10. 

The  centrifugal  motion,  coupled  with  the  volume  expansion  of  the  annular  air  in  the  flame  zone, 
causes  the  data  rate  for  the  annulus  air  seed  to  be  low  enough  that  data  cannot  be  collected  in  a 
region  between  the  hydrogen  jet  and  the  annulus  air  near  the  nozzle  exit.  These  voids  in  the  data 
can  be  noticed  in  Figures  8-9  and  8-10  (up  to  x/D  =  2.65  for  Case  4)  and  in  Figures  8-5  and  8-6 
(up  to  x/D  =  1.06  for  Case  3).  The  radial  velocities  in  the  bulk  of  the  annular  swirling  jet  in  the 
nonreacting  cases  were  less  than  2  m/s  and  directed  inward  [8-6,  8-9]  and  hence  did  not  cause 
such  a  problem  for  LDV  measurements. 

Another  interesting  observation  is  that  although  the  initial  swirl  velocity  is  higher  for  Case  4  than 
Case  3,  the  peak  swirl  velocity  of  the  annular  air  is  approximately  the  same  for  the  two  cases  at 
x/D  =  2.65  as  dictated  by  the  conservation  of  angular  momentum,  due  to  the  larger  radial 
movement  of  the  annular  air  in  the  former  case. 

The  computations  for  Case  4  (Figures  8-9,  8-10,  and  8-11)  are  in  good  agreement  with  data  and 
predict  the  jet  spreading,  the  decay  of  the  swirl,  and  the  temperature  profile  well  for  x/D  =  7.94 
and  beyond.  There  are  some  significant  discrepancies  between  the  data  and  the  computations, 
especially  in  the  mean  temperature  profile  for  x/D  =  2.65.  The  preceding  discussion  of  the 
strong  swirl  and  radial  velocity  and  the  suggestion  by  the  data  in  Figure  8-9  that  the  flow  is 
tending  towards  recirculation  indicate  that  the  boundary  layer  assumptions  are  seriously  violated 
in  the  near-nozzle  region  and,  hence,  the  calculated  radial  and  mean  pressure  gradients  are  in 
error.  Therefore,  it  is  not  reasonable  to  expect  that  the  boundary  layer  algorithm  will  accurately 
calculate  the  flow  in  this  region.  Some  evidence  of  this  is  also  seen  in  Figures  8-5  and  8-6  at  x/D 
=  1.06,  but  the  problem  is  more  severe  for  the  45  deg  swirl  case.  The  large  density  ratio  between 
the  fuel/burnt  gas  and  the  cold  air  also  compounds  the  problem.  It  is  emphasized  that  the 
shortcoming  is  not  that  of  the  pdf  method,  but  rather  it  is  due  to  the  fact  that  the  boundary  layer 
assumptions  are  not  valid  under  the  conditions  encountered  in  this  region.  Elliptic  flow 
calculations  with  the  pdf  method,  to  be  presented  in  the  subsequent  chapters,  is  expected  to  give 
a  better  agreement  in  the  near-nozzle  region.  Indeed,  the  recomputation  of  the  swirling  jet 
hydrogen  flames  using  the  elliptic  pdf  method  presented  in  chapter  10  shows  better  agreement 
with  data  in  the  near-nozzle  region. 


8-14 


10T 


x/D=  26.5 


1.5 

1.0 

.5i 


0+ 


0.0' 


_□ _ o — n  qi  □  D  P  i°  ° — 2_ 


x/D=  15.9 


a  □ 


o  o 


o-oxuxd  n  r  n  D  iQ  .. 


x/D-  7.94 


x/D=  5.29 


x/D=  2.65 


'•°ctr  'i..oTmkoVo  ko  £o  ^  ko 

r/R 


Figure  8-8.  Computed  profiles  of 
Reynolds-averaged  temperature  variance  (lines) 
compared  against  data  (symbols)  for  the  30 
degree  swirl  case  (Case  3). 


Figure  8-9.  Radial  profiles  of  computed 
conditional  mean  axial  velocity  (lines) 
compared  against  data  (symbols)  for  the  45 
degree  swirl  case  (Case  4). 
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Figure  8-10.  Radial  profiles  of  computed 
conditional  mean  tangential  velocity  (lines) 
compared  against  data  (symbols)  for  the  45 
degree  swirl  case  (Case  4). 


Figure  8-11.  Computed  profiles  of 
Reynolds-averaged  mean  temperature  (lines) 
compared  against  data  (symbols)  for  the  45 
degree  swirl  case  (Case  4). 


The  temperature  and  velocity  profiles  for  Case  4  again  indicate  that  the  peak  temperature 
location  falls  at  the  outer  edge  of  the  shear  layer  between  the  jet  and  the  annulus;  however,  it  is 
more  into  the  shear  layer  and  spreads  faster  into  the  shear  layer  than  in  the  two  cases  presented 
previously  (compare  velocity  and  temperature  profiles  at  x/D  =  7.94  for  example). 

Higher  order  turbulence  correlations  are  now  presented  to  compare  the  details  of  the  structure  of 
turbulence  computed  by  the  pdf  method  against  the  data.  Correlations  of  any  order,  including 
velocity-scalar  correlations  (not  presented  here),  can  be  extracted  from  the  pdf  solution  just  as 
velocity  correlations  of  any  order  can  be  extracted  from  the  LDV  data.  Figures  8-12  through 
8-15  present  sample  higher  order  correlations  in  the  developing  region  of  the  flame  for  the  30 
deg  swirl  case  (Case  2).  The  correlations  presented  are  the  turbulent  kinetic  energy  (k),  the 
turbulent  shear  stress  (<uv>),  the  triple  correlation  (<u2v>),  and  the  fourth-order  correlation 
(<v4>),  respectively.  The  correlations  are  in  very  good  agreement  with  the  data  in  terms  of 
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Figure  8-12.  Radial  profiles  of  computed 
conditional  turbulent  kinetic  energy  (lines) 
compare  against  data  (symbols)  for  the  30 
degree  swirl  case  (Case  3). 


Figure  8-13.  Radial  profiles  of  computed 
conditional  turbulent  shear  stress  (lines) 
compared  against  data  (symbols)  for  the  30 
degree  swirl  case  (Case  3). 


magnitudes  and  trends,  and,  as  in  the  nonreacting  case,  show  significant  differences  between  the 
statistics  of  the  fluids  originating  from  different  streams,  especially  between  the  jet  and  the 
annulus  fluids.  These  differences  diminish  as  the  flow  proceeds  downstream  and  the  streams 
mix.  It  should  be  noted  that  the  conventional  second-order  closure  models  compute  only  up  to 
the  second-order  correlations  and  the  third-order  correlations  are  modeled.  Given  that  the  higher 
order  correlations  are  both  difficult  to  measure  and  compute  to  a  high  level  of  accuracy,  the  good 
agreement  between  the  data  and  computations  observed  for  all  quantities  validate  the 
computations  and  the  data. 

The  cases  took  from  approximately  15  minutes  (for  Case  2)  to  23  minutes  (for  Case  4)  of  CPU 
time  on  an  IBM  RS6000/370  workstation.  The  nominal  number  of  particles  used  in  the 
simulations  were  1 10,000  and  the  number  of  spatial  bins  used  (for  averaging)  was  1 10. 
Calculations  with  300,000  particles  yielded  nearly  identical  results.  The  number  of  basis 
functions  used  for  the  spline  fits  was  20. 
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Figure  8-14.  Radial  profiles  of  computed 
conditional  triple  correlation  <u2v>  (lines) 
compared  against  data  (symbols)  for  the  30 
degree  swirl  case  (Case  3). 


Figure  8-15.  Radial  profiles  of  computed 
conditional  fourth  moment  <v4>  (lines) 
compared  against  data  (symbols)  for  the  30 
degree  swirl  case  (Case  3). 


8.7  Summarizing  Remarks 

Computations  using  the  joint  velocity-ffequency-scalar  pdf  method  as  well  as  detailed 
benchmark  quality  measurements  have  been  presented  for  nonswirling  and,  for  the  first  time, 
swirling  hydrogen  jet  diffusion  flames.  The  measurements  and  computations  reported  include 
velocities  (mean  and  higher  moments  up  to  fourth  order)  conditional  upon  the  jet  of  origin  of  the 
fluid  and  temperature  (mean  and  variance)  near  the  burner  exit  and  downstream  locations  up  to 
26.5  jet  diameters.  The  velocities  were  measured  with  a  three-component  LDY  and  the 
temperature  was  measured  using  CARS. 

The  hydrogen  flame,  due  to  the  low  value  of  the  stoichiometric  mixture  fraction  for  hydrogen/air 
combustion,  is  stabilized  at  the  outer  edge  of  the  shear  layer  between  the  hydrogen  and  the 
surrounding  air  jet  as  expected.  As  the  swirl  in  the  air  jet  increases,  the  location  of  the  peak 
temperature  moves  radially  outward  along  with  the  shear  layer  but  still  remaining  at  the  outer 
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edge  of  the  shear  layer.  The  flame  broadens  and  spreads  more  into  the  shear  layer  due  to 
enhanced  mixing  due  to  swirl.  Other  relative  features  of  the  flames  have  been  studied  and 
discussed  in  detail. 

The  computations  are  in  good  agreement  with  data  including  those  for  fourth-order  turbulent 
correlations,  demonstrating  and  further  validating  the  ability  of  the  pdf  method  to  accurately 
compute  the  details  of  the  flow  and  the  local  turbulence  structure.  An  important  finding  is  that 
although  Favre  (density-weighted)  averages  are  recommended  and  typically  used  in  the 
computation  of  reacting  flows,  the  conventional  Reynolds  average  is  the  appropriate  average  to 
compare  against  temperature  data  from  CARS. 

A  relatively  simple  boundary-layer  algorithm  with  a  relatively  simple  model  for  thermochemistry 
used  for  the  computations  provided  a  remarkably  good  agreement  with  velocity  and  temperature 
data,  except  very  close  to  the  nozzle  for  the  high  swirl  case  where  the  boundary-layer 
assumptions  are  not  likely  to  be  valid.  Given  the  current  agreement  with  the  data,  further 
complexities  in  the  modeling  of  the  thermochemistry  such  as  multistep  chemistry,  differential 
diffusion  of  species,  and  nonunity  Lewis  number  effects  may  not  be  warranted  for  the  flames 
studied.  The  computations  took  a  maximum  of  23  minutes  on  an  IBM  RS6000/370  workstation. 
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9.  BACKWARD  FACING  STEP  FLOW 


9.1  Overview 


The  algorithms  to  determine  the  mean  pressure  field  for  elliptic  flow  calculations  with  the  pdf 
method  (sections  5.2  and  5.3)  are  applied  and  validated.  Previous  computations  of  elliptic  flows 
with  the  pdf  method  have  involved  accompanying  conventional  finite-volume  based  calculations 
to  provide  the  mean  pressure  field.  The  algorithms  developed  and  here  permit  the  mean  pressure 
field  to  be  determined  within  the  pdf  calculations.  The  pdf  method  incorporating  the  pressure 
algorithm  is  applied  to  the  flow  over  a  backward-facing  step.  The  results  are  in  very  good 
agreement  with  data  for  the  reattachment  length,  mean  velocities,  and  turbulence  quantities 
including  triple  correlations. 


9.2  Background 

As  discussed  in  section  2.3.  there  are  no  theoretical  limitations  for  applying  the  pdf  method  for 
complex  (2-D  or  3-D)  elliptic  flows.  However,  a  suitable  algorithm  for  extracting  the  mean 
pressure  field  within  the  pdf  solution  algorithm  is  needed.  A  pdf  method  for  elliptic  flows  was 
first  demonstrated  by  Anand  et  al.  [9-1].  In  this  method,  the  Monte  Carlo  (MC)  calculations  for 
the  pdf  are  combined  with  conventional  finite-volume  (FV)  calculations  for  mean  fields  (Ref. 
11).  The  mean  pressure  field  and  the  turbulent  time  scale  are  supplied  to  the  MC  calculations  by 
the  FV  calculations.  In  turn,  the  MC  calculations  supply  to  the  FV  calculations  the  Reynolds 
stresses  extracted  from  the  pdf  solution,  thereby  avoiding  the  use  of  conventional  turbulence 
models.  This  combined  algorithm  exploits  the  advantages  of  both  methods  to  yield  an 
economical  algorithm  for  pdf  calculations  of  complex  elliptic  flows.  This  method  was  intended 
to  demonstrate  the  feasibility  of  pdf  calculations  for  such  flows. 

However,  for  more  complex  flows,  especially  variable  density  and  reacting  flows  with  fast  or 
finite-rate  multi-step  chemistry,  the  coupling  between  the  two  methods  becomes  quite  complex 
and  it  becomes  more  difficult  to  fully  exploit  the  advantages  of  the  pdf  method.  Therefore,  it  is 
desirable  to  perform  pdf  calculations  independently,  without  recourse  to  FV  calculations.  This 
chapter  describes  the  application  of  two  algorithms  developed  through  which  pdf  calculations  are 
performed  independently  and  economically  (in  terms  of  computer  resources  required)  for 
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complex  elliptic  flows.  In  addition,  the  treatment  of  wall  boundaries  is  included  in  the 
calculations.  The  method  is  applied  to  a  sample  elliptic  flow,  namely  the  recirculating  flow 
behind  a  backward-facing  step  and  compared  against  data  of  Pronchick  and  Kline  [9-2]. 


9.3  Test  Case  and  Initial  Conditions 

The  schematic  of  the  flow  is  shown  in  Figure  9-1.  The  step  height  is  H  and  the  centerline  mean 
axial  velocity  in  the  channel  is  Uref.  The  fluid  used  in  the  experiments  was  water,  and  the 
Reynolds  number  based  on  H  and  Uref  was  approximately  12,000. 

The  inlet  velocity  pdf  was  prescribed  to  be  joint-normal  with  the  means  and  covariances  of 
the  velocity  components  being  those  of  fully  developed  turbulent  flow  in  the  inlet  channel  as 
calculated  by  the  k-e  model.  The  inlet  turbulent  kinetic  energy  k  is  distributed  into  the  three 
components  according  to  k  =  <u2>  =  2<v2>  =  2<w2>.  The  inlet  shear  stress  <uv>  is  prescribed 
to  be  linear  with  a  value  of  zero  at  the  centerline  of  the  inlet  channel  and  a  peak  value  of  0.03k 
near  the  wall  of  the  inlet  channel.  The  turbulence  levels  at  the  inlet  have  little  effect  on  the  final 
results  since  most  of  the  shear  and  turbulence  production  occurs  in  the  recirculation  zone  within 
the  solution  domain.  The  initial  pdf  in  the  solution  domain  was  set  equal  to  the  inlet  pdf  and  the 
initial  mean  pressure  was  zero  everywhere. 

Computations  were  performed  using  both  the  elliptic  flow  pressure  algorithms  described  in 
Chapter  5.  The  latest  models  for  velocity  and  frequency  available  at  the  time  that  these 
algorithms  were  developed,  were  used.  Results  are  presented  in  the  next  two  subsections. 


*R 


H  -  0.0762  m  W3/(W3-H)«  1.43 

W3  -  0.254  m  Uref  -  0.196  m/s 


Figure  9-1.  Schematic  of  the  flow  considered.  Backward-facing  step  studied  experimentally 

by  Pronchick  and  Kline  (1983). 
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9.4  Computations  with  Pressure  Algorithm  I 


9.4.1  Computational  Details 

The  pressure  algorithm  described  briefly  in  section  5.2  is  applied  and  validated.  There  are 
considerable  details  in  the  formulation  and  implementation  of  the  algorithm.  For  the  sake  of 
brevity  they  are  not  presented  here.  The  reader  if  referred  to  Anand  et  al.  [9-3]  for  details. 

The  simple  Langevin  model  (SLM  I)  is  used  for  the  velocities.  The  stochastic  turbulence 
frequency  models  were  still  being  developed  at  the  time  the  PA  I  algorithm  was  developed  and 
applied.  Hence  the  mean  turbulence  frequency  was  supplied  to  the  pdf  calculations  in  the 
following  manner.  The  turbulent  frequency  (<co>  =  <£>/k,  where  <e>  is  the  mean  dissipation 
rate  of  k)  is  calculated  at  each  step  from  the  current  value  of  k  by 

<co>  =  k0.5//  ,  (9.1) 

where  the  turbulent  length  scale  /  is  supplied  from  the  k-e  solution  of  Ref.  [9-1]  (see  Eq.  19),  and 
is  held  constant  throughout  the  calculations: 

/  =  kL5/<e>  .  (9.2) 

The  choice  of  supplying  /  was  preferred  over  the  choices  of  supplying  <(D>  or  <e>  from  the  k-e 
solution  since  the  profiles  of  /  were  nearly  uniform  across  the  flow  (approximately  0.02  13) 
except  near  the  walls  where  the  values  decrease.  This  way,  the  /  used  in  the  calculations  (Eq. 
9.1)  would  be  more  consistent  with  the  current  values  of  k. 

The  number  of  particles  used  in  the  simulation  was  60,000.  The  computations  were  performed 
for  250  steps  with  pressure  updates  every  20  steps  and  the  pressure  relaxation  factor  of  0.8.  The 
time  interval  for  each  step  was  0.3  H/Uref  (see  Figure  9-1  for  the  definition  of  H  and  Uref).  This 
time  interval  was  less  than  one-tenth  of  the  turbulent  time  scale  over  most  of  the  flow  domain 
except  for  small  regions  near  the  wall  where  it  was  between  one-fifth  and  one-tenth  of  the 
turbulent  time  scale.  The  calculations  reached  a  steady-state  around  the  200th  step.  The  final 
results  shown  here  (except  for  mean  pressure  which  is  taken  from  the  last  step)  are  calculated 
from  spline-fits  to  sums  formed  over  50  steps  after  steady-state  has  been  reached.  This  has  an 
effect  similar  to  having  3  million  (60,000  x  50)  particles  in  the  simulation  thereby  minimizing 
the  statistical  error  in  the  results  which  is  inversely  proportional  to  the  square  root  of  the  number 
of  particles. 
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9.4.2  Results  and  Discussion 


The  computed  reattachment  length  of  6.5H  is  in  good  agreement  with  the  experimentally 
observed  reattachment  length  (xre)  of  6.75H.  It  should  be  noted  that  this  agreement  is  better 
than  achieved  by  any  finite-volume  calculations. 

The  results  for  mean  velocities  and  various  turbulent  quantities  are  shown  in  Figure  9-2  through 
9-10.  The  results  are  shown  for  different  axial  stations  normalized  with  respect  to  the 
experimental  reattachment  length: 

X*  =  (X  -  xre)  /xre  .  (9-3) 

Overall,  the  results  are  in  excellent  agreement  with  data  both  in  magnitude  and  qualitative  trends. 
It  is  noteworthy  that  the  turbulence  quantities,  especially  the  third  moments  (which  cannot  be 
calculated  with  the  k-e  model,  and  are  modeled  in  second-moment  closures)  are  in  good 
agreement  with  data.  It  is  also  to  be  noted  that  the  turbulence  quantities  are  an  order  of 
magnitude  lower  than  the  reference  velocity  Uref. 

The  contours  of  the  calculated  mean  pressures  are  shown  in  Figure  9-11.  The  magnitudes  and 
shapes  of  the  contours  are  comparable  to  those  obtained  from  previous  calculations  with  the  k-e 
model  and  combined  FV/MC  pdf  method  [9-1].  This  fact,  in  conjunction  with  the  accurate 
prediction  of  the  mean  velocity  field,  shows  that  the  pressure  algorithm  is  working  satisfactorily. 

Approximately  1.3  Mega  words  of  storage  was  required  for  the  computations.  The  total  CPU 
time  required  for  the  computations  (250  steps)  on  a  CRAY  XMP-22  was  approximately  11 
minutes  (about  6  minutes  on  an  IBM  RS6000/370  and  about  3  minutes  on  a  CRAY  C90).  This 
reflects  an  order  of  magnitude  reduction  over  the  estimated  time  it  would  have  taken  with  the 
initial  straight  forward  implementation  of  the  algorithm.  Nearly  55  percent  of  the  computational 
time  is  spent  in  computations  related  to  the  pressure  algorithm,  namely  forming  sums  and  splines 
for  pressure  updates  and  velocity  corrections,  in  spite  of  limiting  the  frequency  of  pressure 
updates  and  improving  the  efficiency  of  the  spline-fitting  routines.  However,  the  total 
computational  time  of  approximately  1 1  minutes  and  the  storage  required  are  still  very  modest 
for  pdf  calculations  of  an  elliptic  flow. 
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Figure  9-10.  Computed  <v2u>  profiles  (lines)  compared  against  data  (symbols) 
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Figure  9-11.  Computed  contours  of  mean  pressure  (normalized  by  pUref  ) 


9.4.3  Summarizing  Remarks 


An  algorithm  to  determine  the  mean  pressure  in  steady-state  constant  or  variable-density  elliptic 
flow  calculations  with  the  pdf  method  has  been  developed  and  implemented  economically.  The 
algorithm  and  the  pdf  calculations  have  been  tested  for  a  sample  recirculating  flow;  namely  the 
flow  behind  a  backward-facing  step. 

The  computed  results  are  in  excellent  agreement  with  data  for  the  mean  velocity  field  and 
various  turbulence  quantities  including  triple  corrections.  The  computed  reattachment  length  is 
in  excellent  agreement  with  the  experimental  value.  These  comparisons  have  served  to  validate 
the  pressure  algorithm  and  its  implementation. 

9.5  Computations  with  Pressure  Algorithm  II 

The  backward  facing-step  flow  was  also  computed  using  the  newer  pressure  algorithm  (PA  II, 
section  5.3).  As  explained  in  that  section,  the  newer  algorith  is  more  robust  and  is  easier  to 
extend  to  complex  geometry  and  body-fitted  grids.  The  models  used  included  the  latest  fequency 
model  (SFM  II)  and  the  corresponding  velocity  model  (SLM  II).  The  results  (not  shown  here) 
were  nearly  the  same  as  the  results  presented  in  the  previous  section  (Figures  9-2  through  9-11) 
and  the  were  in  excellent  agreement  with  data.  The  computed  reattachment  length  of  6.77H  was 
in  excellent  agreement  with  the  experimentally  observed  reattachment  length  (x^g)  of  6.75H. 

These  results  validate  the  new  pressure  algorithm  and  the  models  for  a  constant-density  flow. 
The  validation  for  swirling,  recirculating  and  reacting  flows  is  presented  in  the  next  chapter. 
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10.  SWIRLING  AND  RECIRCULATING  FLOW  COMBUSTORS 


10.1  Overview 


Calculations  are  reported  for  recirculating  swirling  reacting  flows  using  the  joint  velocity- 
frequency-scalar  pdf  method.  The  pdf  calculations  reported  here  are  based  on  a  newly  developed 
solution  algorithm  for  elliptic  flows  (PA  II),  and  on  newly  developed  models  for  turbulent 
frequency  and  velocity  (SFM  II  and  SLM  II)  that  are  simpler  than  the  SFM  I  and  RLM  models 
used  in  some  of  the  calculations  reported  in  the  previous  chapters.  Calculations  are  performed 
for  two  different  gas-turbine-like  swirl  combustor  flows  for  which  detailed  measurements  were 
made  at  WPAFB.  The  computed  results  are  in  good  agreement  with  experimental  data.  More 
details  of  the  study  are  presented  in  Ref.  10-1. 


10.2  Background 

The  present  study  represents  the  first  fully  self-contained  pdf  calculations  for  elliptic  reacting 
flows  and  incorporates  the  elliptic  flow  solution  algorithm  as  well  as  the  stochastic  frequency 
model.  With  a  view  to  making  the  method  more  robust,  easier  to  implement  and  affordable  for 
complex  multidimensional  flows,  a  significantly  different  elliptic  flow  algorithm  (or  pressure 
algorithm)  than  PA  I  has  been  developed  and  implemented.  The  models  for  turbulent  frequency, 
scalar  mixing,  and  velocity  have  also  been  considerably  simplified.  (See  chapters  3, 4  and  5.) 

The  pdf  method  is  validated  against  benchmark  experimental  data  for  two  gas-turbine-like  swirl 
combustors  and  also  compared  against  previous  pdf  solutions  for  one  of  the  combustors 
considered. 


10.3  Computational  Details 

10.3.1  The  PDF  Method  --  Modeling  and  Solution  Algorithm 

As  mentioned  above,  the  PA  II  algorithm  for  pressure,  the  SFM  II  model  for  frequency,  the  SLM 
II  model  for  velocity,  and  the  IEM  model  for  scalars  are  used  in  the  study  presented  here.  One  of 
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of  the  features  of  the  new  algorithm  is  the  use  of  a  different  (less  computer  intensive  but  more 
sophisticated)  technique  for  estimating  means  and  smoothing  the  profiles  then  the  regular 
ensemble  averaging  and  cubic  B-splines  used  in  the  boundary-layer  and  PA  I  algorithms.  In  the 
present  study,  means  of  any  functions  of  the  independent  variables  are  determined  by  a 
sophisticated  ensemble  averaging  procedure  (cloud-in-cell  estimate  using  bi-linear  basis 
functions)  followed  by  smoothing  using  local  linear  least  squares  [10-2].  Additionally  time¬ 
averaging,  with  a  low  time  constant  initially  and  a  higher  value  near  convergence,  is  used  for 
mean  quantities  to  further  reduce  the  statistical  error.  The  details  of  the  pressure  algorithm  can 
be  found  in  Chapter  5.3  and  in  Ref.  [10-3]. 

The  time  increment.  At,  for  each  step  is  chosen  to  be  a  fraction  (=0.1)  of  the  minimum  of  a)  the 
inverse  of  the  maximum  mean  turbulence  frequency  in  the  computational  domain  or  b)  the 
minimum  characteristic  time  for  any  particle  to  cross  a  computational  cell  based  on  the  mean  and 
variance  of  velocity  in  the  cell.  All  the  particles  in  the  computation  are  marched  with  the  same 
time  increment.  The  particle  evolution  equations  are  integrated  over  the  time  step  with  an 
accuracy  of  second-order  or  better. 

It  should  be  noted  that  the  models  described  above  for  velocity,  frequency  and  scalar  mixing  are 
all  being  used  in  the  joint  pdf  method  for  the  first  time  to  compute  general  (inhomogeneous, 
swirling,  recirculating)  turbulent  reacting  flows.  As  such  the  present  study  serves  to  validate  the 
elliptic  flow  algorithm  as  well  as  the  models  used. 


10.3.2  Thermochemistry 

Hydrogen  and  methane  flames  are  studied  in  the  present  work.  A  fast  equilibrium  chemistry 
model  is  used  for  the  hydrogen  flame  calculations  because  the  time  scale  for  hydrogen-air 
reaction  is  very  small  compared  to  the  turbulent  time  scale.  For  the  hydrogen  case,  the  only 
scalar  variable  in  the  calculations  is  the  mixture  fraction.  Temperature  is  also  included  but  is 
needed  for  output  only.  The  mixture  fraction  is  a  conserved  variable  (reaction  rate  is  zero).  The 
density  and  temperature  are  determined  as  equilibrium  properties  from  the  mixture  fraction. 

For  methane  flame  calculations,  a  general  2-step  chemistry  by  Westbrook  and  Dryer  [10-4]  for 
saturated  hydrocarbon  fuels  is  used.  The  two  steps  are: 


(10.1) 
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(10.2) 


co+|o2  =  co2 . 

In  addition  to  mixture  fraction,  two  more  scalar  variables— namely  the  mass  fractions  of  carbon 
dioxide  and  water— are  included  in  the  pdf  calculations.  The  temperature  and  density  are 
determined  as  functions  of  these  three  scalar  variables. 

For  both  the  fast  chemistry  and  the  2-step  chemistry  models,  lookup  tables  were  created  to 
reduce  the  CPU  requirements  of  the  calculations.  In  the  case  of  the  fast  chemistry,  a  one¬ 
dimensional  table  is  created,  and  for  the  2-step  chemistry,  a  three-dimensional  table  is  generated. 
For  the  2-step  chemistry  calculations,  the  table  is  generated  for  a  given  specific  time  increment. 
At,  used  by  the  flow  calculations  (At  =  2.5x1 0'^  s  in  the  present  calculations).  In  the  table 
generation  processes,  the  NASA  CEC  thermal  data  were  used  to  calculate  the  variable  specific 
heats  and  the  temperature. 


10.4  Results  and  Discussion 


Results  are  presented  for  two  laboratory  swirl  combustor  configurations  which  have  the  essential 
flow  features  of  gas  turbine  combustors,  namely  swirl,  recirculation,  large  velocity  gradients, 
turbulence  and  combustion.  The  experiments  were  conducted  by  researchers  at  the  University  of 
Dayton  Research  Institute  at  Wright  Patterson  Air  Force  Base,  Ohio.  The  velocity  measurements 
were  made  using  a  three-component  laser  Doppler  velocimeter  (LDV)  and  temperature 
measurements  were  made  using  coherent  anti-Stokes  Raman  spectroscopy  (CARS).  (See  chapter 
8  for  more  details.) 


10.4.1  Swirling  Hydrogen  Diffusion  Flame 

Figure  10-1  shows  the  schematic  of  the  swirling  jet  diffusion  flame  combustor  configuration. 
The  test  case  considered  had  a  central  fuel  (hydrogen)  jet  bulk  velocity  of  100  m/s,  the  swirling 
air  bulk  velocity  of  20  m/s  and  the  nonswirling  coflow  air  velocity  of  4  m/s.  The  swirler  vane 
angle  was  30  deg.  and  swirl  number  for  the  swirling  jet,  calculated  from  the  measured  velocities, 
was  0.382.  Detailed  measurements  for  mean  velocities  and  turbulent  correlations  up  to  fourth 
order  are  reported  (see  Ref.  10-5  and  references  listed  therein)  at  several  downstream  locations 
starting  from  the  axial  location  x  =1.5  mm  from  the  nozzle.  While  no  species  measurements 
were  made,  mean  and  variance  of  the  temperature  are  reported  at  the  same  locations. 
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AIR 

Figure  10-1.  Schematic  of  the  hydrogen  swirling  jet  diffusion  flame  configuration. 

The  flow  was  previously  calculated  by  Anand  et  al.  [10-5]  using  a  joint  pdf  method.  (The  study 
is  reported  in  Chapter  8.)  Since  the  flow  is  primarily  parabolic  in  nature  (with  no  recirculation), 
the  pdf  solution  algorithm  was  based  on  boundary-layer  assumptions  with  extensions  for  swirling 
flows.  The  method  also  used  more  sophisticated  models,  namely  the  stochastic  frequency  model 
SFM  I  in  conjunction  with  the  refined  Langevin  model  RLM  for  velocity  in  which  the 
instantaneous  particle  frequency  rather  than  a  mean  frequency  is  used  in  the  random  term  (see 
sections  3.3  and  3.4).  The  reason  for  calculating  the  swirling  hydrogen  flames  again  with  the 
present  method  is  to  not  only  validate  the  method  and  the  models  but  also  to  assess  if  the  elliptic 
flow  algorithm  can  better  resolve  the  flow  in  regions  where  boundary-layer  assumptions  (e.g., 
neglect  of  axial  gradients)  are  questionable  (see  discussion  in  section  8.5). 

The  present  computations  were  performed,  on  an  IBM  RS6000/370,  using  a  nonuniform  grid 
(31  in  x  by  61  along  the  radius  r)  with  about  190  particles  per  cell.  Increasing  the  nominal 
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number  of  particles  per  cell  to  290  produced  nearly  the  same  results.  The  inlet  boundary 
conditions  were  taken  from  experimental  data  as  described  in  Ref.  10-5.  Calculations  were 
performed  for  2000  time  steps. 

In  the  figures  to  be  presented  for  this  case  the  measured  mean  axial  velocity  on  the  centerline  at 
the  nozzle  exit,  <U>0C,  which  is  130.3  m/s  in  this  case,  is  used  to  normalize  the  velocity 
statistics.  The  axial  distances  and  radius  are  normalized  by  the  nozzle  diameter,  D  (=  9.45  mm) 
and  nozzle  radius  R  (=  D/2)  respectively.  The  temperature  results  are  normalized  by  the 
stoichiometric  temperature,  Ts^  (=  2377  K).  For  the  experimental  data  presented,  the  open 
squares  represent  data  conditioned  on  the  inner  fuel  jet,  the  solid  circles  represent  data 
conditioned  on  the  swirling  air  jet,  and  the  inverted  triangles  represent  the  data  conditioned  on 
the  outer  coflow  air. 

Figure  10-2  shows  the  convergence  history  for  the  (unnormalized)  mean  axial  and  radial 
velocities  at  the  monitoring  location  x/D  =  10  and  r/R  =  4.  The  figure  shows  that  the  solution  has 
converged  and  a  steady-state  has  been  reached.  Typically  oscillations  are  seen  during  the  initial 
(first  few  hundred)  steps  before  the  solution  settles  down  and  reaches  a  steady-state.  Also, 
steady- state  and  convergence  were  usually  achieved  by  about  1500  steps  for  both  the  cases 
presented  here. 


Figure  10-2.  Convergence  history  for  the  mean  axial  and  radial  velocities  for  the  hydrogen 

flame  computations. 
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Figure  10-3  shows  the  radial  profiles  of  normalized  mean  axial  velocity  <U>  at  different 
downstream  locations  compared  against  data.  The  measurements  are  conditional  on  the  origin  of 
the  fluid  and  are  made  by  seeding  (for  LDV)  each  of  the  jets  (fuel,  annulus  and  coflow) 
individually.  Differences  in  the  velocity  statistics  for  each  of  the  jets  can  be  seen  in  the  data. 
The  pdf  method  can  calculate  these  conditional  values  without  requiring  additional  modeling, 
and  such  calculations  have  been  presented  in  previous  chapters.  However,  in  the  present  study, 
only  the  unconditional  quantities  are  calculated  and  compared  against  data.  The  present  results 
(solid  lines)  are  also  compared  against  results  from  the  boundary-layer  (b-1)  algorithm  (dashed 
line).  Figure  10-3  shows  that  the  present  results  are  in  excellent  agreement  with  data  at  all 
stations.  Also,  the  present  calculations  are  in  better  agreement  with  data  than  the  b-1  results  not 
only  with  respect  to  the  overall  spreading  but  also  in  the  near-nozzle  region  (x/D  =  1.06  and 
2.65)  where  the  flow  is  tending  toward  recirculation  (near  r/R  =  1.5)  and  the  b-1  approximations 
are  inaccurate. 

Similar  observations  can  be  made  for  the  swirl  velocity  results  presented  in  Figure  10-4. 
Although  the  b-1  calculations  are  in  good  agreement  with  data,  the  present  calculations  show  a 
better  agreement. 

The  profiles  of  mean  (Reynolds-averaged)  temperature,  T,  presented  in  Figure  10-5  show  that 
the  transport  and  mixing  of  the  fuel  is  well  calculated  in  the  present  study,  resulting  in  a  very 
good  agreement  with  the  temperature  data.  (Note  that  CARS  measurements  are  closer  to 
Reynolds-averaged  values  than  density-weighted  values  as  discussed  in  section  8.5.)  The  present 
results  are  better  than  the  b-1  results  at  the  downstream  locations,  but  near  the  nozzle  (e.g.,  x/D  = 
2.65)  the  present  results  show  a  lower  peak  and  a  greater  spread  than  the  data  and  the  b-1  results 
show.  Results  for  the  temperature  variance  from  both  the  methods  (not  shown  here)  were  overall 
in  good  agreement  with  data  although  some  differences  consistent  with  their  mean  temperature 
profiles  were  observed. 

Profiles  for  turbulent  kinetic  energy  presented  in  Figure  10-6  show  that  both  the  calculations  are 
in  good  agreement  with  data  in  the  region  up  to  x/D  =  5.29  while  the  b-1  calculations  overpredict 
the  kinetic  energy  at  downstream  locations  which  is  consistent  with  the  lack  of  spreading  in  the 
mean  velocity  profiles.  Sample  results  from  the  present  calculations  for  third  and  fourth  order 
turbulent  correlations  presented  in  Figure  10-7  are  in  good  agreement  with  data. 

Overall,  the  results  are  in  very  good  agreement  with  data,  and  are  as  good  or  better  than  those 
obtained  with  the  boundary-layer  calculations. 
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Figure  10-3:  Computed  mean  axial  velocity  Figure  10-4:  Computed  mean  swirl  velocity 
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Figure  10-5:  Computed  mean  temperature 
profiles  compared  against  data  for  the 
hydrogen  flame. 


Figure  10-6:  Computed  turbulent  kinetic 
energy  profiles  compared  against  data  for 
the  hydrogen  flame.. 


3] 


Figure  10-7:  Computed  profiles  of  higher-order  turbulent  correlations  compared  against  data  for 

the  hydrogen  flame. 
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10.4.2  Methane  Step-Swirl  Combustor 


The  step-swirl  combustor  shown  in  Figure  10-8  is  an  extension  of  the  jet  diffusion  flame 
combustor  (Figure  10-1),  and  is  closer  to  a  practical  gas  turbine  combustor.  It  consists  of  a 
central  air  jet  (20  mm  diameter)  surrounded  by  an  annular  fuel  tube  (29  mm  o.  d.)  which  is  again 
surrounded  by  a  swirling  air  jet  (40  mm  o.  d.  taken  to  be  the  characteristic  diameter  D). 
Measurements  for  this  case  are  reported  in  Ref.  [10-6].  For  the  case  considered  here  the  inner  air 
jet  was  nonswirling  and  the  outer  air  had  a  30  deg.  vane-angle  swirler.  The  bulk  velocity  for  the 
inner  air,  fuel  and  outer  air  jets  are  14.4, 2.5,  and  8.6  m/s,  respectively. 

The  velocity  data  reported  for  this  combustor  are  also  conditional  velocities.  Unfortunately,  the 
authors  [10-6]  were  unable  to  measure  the  velocities  conditional  on  the  outer  air  jet  due  to 
practical  difficulties  such  as  the  LDV  seed  particles  striking  the  optical  windows  and  clogging 
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Figure  10-8.  Schematic  of  the  step-swirl  combustor.  The  fuel  used  is  methane. 
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them  up.  The  outer  swirling  flow  has  a  major  effect  on  the  development  of  the  flow  and  it  is 
crucial  to  have  accurate  inlet  conditions  to  accurately  simulate  the  flow.  The  computations  also 
show  high  sensitivity  to  the  inlet  profiles,  especially  since  the  comparison  with  data  is  made  in 
the  region  x/D<2  up  to  which  measurements  were  made.  In  the  present  study  the  inflow  velocity 
profiles  had  to  be  reconstructed  based  on  existing  experimental  data  and  the  overall  mass  flow 
rates  through  the  different  streams. 

The  computations  were  performed  using  a  nonuniform  41  x  41  grid  using  200  particles  per  cell. 
The  computations  were  made  for  2000  steps  with  convergence  achieved  in  most  of  the  flow  field 
by  about  1500  steps.  The  results  are  presented  in  Figures  10-9  through  10-14.  The  values  used 
for  normalization  in  the  figures  are  <U>0C  =  21.6  m/s,  Tst  =  2272  K,  D  =  40  mm,  and  R  =  D/2. 
For  the  experimental  data  presented,  the  solid  symbols  represent  data  conditioned  on  the  inner  air 
jet  and  the  open  symbols  represent  data  conditioned  on  the  fuel  jet. 

The  profiles  of  mean  axial  velocity,  <U>,  presented  in  Figure  10-9  show  the  calculations  capture 
the  overall  flow  features  well.  Although  the  recirculation  is  well  predicted,  the  location  and 
radial  extent  of  the  recirculation  zone  which  are  very  sensitive  to  the  inlet  mean  radial  and  swirl 
velocities  assumed  for  the  outer  swirling  jet,  are  underpredicted. 

The  profiles  of  mean  radial  velocity,  <V>,  in  Figure  10-10  show  the  expected  trends  although 
data  are  not  available  in  the  critical  regions  where  the  largest  radial  velocities  are  present  Note 
that  computed  results  are  unconditional  and  are  expected  to  be  lower  than  the  fuel  conditioned 
radial  velocity  at  the  outer  edge  of  the  fuel  jet  as  seen  at  x/D  =  0.5.  Figure  10-11  shows  that  the 
mean  swirl  (or  tangential)  velocity,  <W>,  is  well  predicted  both  in  terms  of  the  peak  location  and 
the  decay  downstream. 

The  mean  temperature  profiles  presented  in  Figure  10-12  predict  the  shapes  of  the  measured 
profiles  well  and  for  the  most  part  agree  in  magnitude  with  the  data.  The  profile  of  the  fuel  mass 
fraction  at  the  inlet  significantly  influences  the  temperature  distribution  at  the  near-nozzle 
locations  at  which  comparisons  are  being  made.  Although  the  fuel  tube  only  supplies  fuel, 
considerable  mixing  takes  place  even  as  the  fuel  is  leaving  the  fuel  tube,  and  an  assumption  of  a 
plug  flow  profile  leads  to  a  much  worse  comparison  with  temperature  than  with  a  smooth  but 
sharply  peaked  profile  assumed  for  the  computations  shown. 

The  profiles  of  turbulent  kinetic  energy  and  a  fourth  order  turbulent  correlation  shown  in  Figures 
10-13  and  10-14,  respectively,  are  in  reasonably  good  agreement  with  data.  Overall,  the  results 
are  in  good  agreement  with  data  for  all  the  quantities  considering  the  uncertainty  in  the  inlet 
conditions. 
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Figure  10-9:  Computed  mean  axial  velocity 
profiles  (lines)  compared  against  data 
(symbols)  for  the  methane  flame. 


Figure  10-10:  Computed  mean  radial 
velocity  profiles  (lines)  compared  against 
data  (symbols)  for  the  methane  flame. 
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Figure  10-13:  Computed  turbulent  kinetic 
energy  profiles  (lines)  compared  against 
data  (symbols)  for  the  methane  flame. 
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Figure  10-14:  Computed  profiles  of  fourth 
moment  of  axial  velocity  (lines)  compared 
against  data  (symbols)  for  the  methane 
flame 
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The  results  for  the  hydrogen  and  methane  cases  have  validated  the  new  models  and  the  elliptic 
flow  algorithm  used.  The  calculations  represent  the  first  quantitative  results  from  the  new  code 
incorporating  the  algorithm  and  models.  The  results  compare  very  well  with  the  detailed  data 
from  practical  combustors. 


10.5  Summarizing  Remarks 

Computations  using  the  joint  pdf  approach  have  been  reported  for  two  swirl  combustor 
configurations.  The  study  uses  a  newly  developed  solution  algorithm  for  elliptic  flows  and  new 
simplified  models  for  velocity  and  turbulence  frequency.  The  methane  combustor  calculations 
represent  the  first  fully  self-contained  joint  pdf  calculations  for  elliptic  reacting  flows.  The 
results  for  both  combustors  are  in  good  agreement  with  data.  The  study  serves  to  further  validate 
the  joint  pdf  method  and  the  models,  and  is  a  significant  step  in  the  development  of  a  pdf-based 
combustor  design  system. 

The  ability  of  the  joint  pdf  method  to  accurately  calculate  the  mean  and  turbulent  velocity  fields, 
scalar  transport,  and  temperature  using  multi-step  finite-rate  chemistry  offers  significant 
advantages  for  its  use  in  the  design  of  current  and  future  high  performance  and  low  emissions 
gas  turbine  combustors. 

The  present  results  are  compared  against  calculations  using  the  scalar  pdf  method  (in  which  the 
joint  pdf  of  only  the  scalars  is  considered)  and  other  conventional  turbulent  combustion  models 
in  the  next  chapter.  The  study  demonstrates  the  advantages  and  the  superior  accuracy  of  pdf 
methods,  in  particular  the  joint  velocity-scalar  pdf  method. 
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11.  ASSESSMENT  OF  PDF  VERSUS  CONVENTIONAL  MODELS 


11.1  Overview 


The  objective  of  the  work  presented  in  this  chapter  is  to  systematically  assess  the  performance  of 
the  joint  pdf  method  against  several  other  methods  and  models  used  in  practical  combustion 
computations.  The  methods  evaluated  include  laminar  chemistry,  the  eddy  breakup  model,  the 
assumed  pdf  model,  the  scalar  pdf  method,  and  the  joint  velocity-scalar  pdf  method.  As 
explained  earlier,  the  treatment  of  turbulent  transport  and  the  two-way  turbulence/chemistry 
interactions  are  some  of  the  key  issues  in  the  computations  of  turbulent  reacting  flows.  The 
assessment  is  done  by  comparing  the  different  model  results  against  benchmark  data  from  two 
different  gas-turbine-like  combustors  presented  in  the  previous  chapter.  The  present  work  shows 
that  the  pdf  methods,  especially  the  joint  pdf  method,  offer  considerable  advantage  in  dealing 
with  these  key  issues,  and,  as  a  result,  are  more  accurate  than  the  conventional  models.  The 
study  presented  in  this  chapter  was  conducted  under  Allison’s  IR&D  program.  It  is  reported  in 
detail  in  Ref.  11-1.  An  application  of  the  pdf  method  to  a  practical  design  problem  and  the 
advantages  offered  are  the  also  presented  and  discussed. 


11.2  Background 

As  discussed  in  Chapter  2,  important  physical  phenomena  to  be  considered  in  the  simulation  of 
turbulent  reacting  flows  include  the  chemistry  closure  problem,  turbulent  transport  and  the 
two-way  interaction  between  turbulence  and  chemistry.  The  chemistry  closure  problem  is  caused 
by  the  following  two  factors:  (1)  according  to  the  Arrhenius  law,  the  reaction  rate  coefficient  is 
an  exponential  function  of  the  temperature  and  (2)  in  turbulent  flames,  large  fluctuations  in 
temperature  and  species  concentrations  occur.  Because  of  these  two  factors,  the  mean  reaction 
rate  is  often  very  different  from  the  reaction  rate  calculated  from  the  mean  temperature  and 
species  concentrations,  and  the  modeling  of  this  term  becomes  a  major  hurdle  in  turbulent 
combustion  simulations.  This  effect  of  turbulence  on  chemistry  is  part  of  the  two  way 
turbulence/chemistry  interaction.  In  the  other  direction  of  this  interaction,  the  fluctuation  in 
reaction  rates  causes  density  fluctuations  in  the  turbulent  flames,  which  in  turn  causes  turbulence 
modulation  and  other  effects  including  counter-gradient  diffusion. 
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The  pdf  methods,  including  the  scalar  pdf  and  joint  velocity-scalar  pdf  methods,  offer  key 
advantages  for  addressing  these  issues.  In  the  present  study,  we  try  to  assess  the  performances, 
including  the  accuracy  and  computer  resources  requirements,  of  these  methods,  and  compare 
them  with  the  conventional  combustion  models  such  as  the  assumed  pdf  model,  the  eddy  breakup 
model,  and  laminar  chemistry.  The  accuracies  of  these  methods  are  assessed  by  comparing  the 
numerical  results  for  a  swirling  hydrogen  jet  flame  and  a  methane  step-swirl  combustor  against 
benchmark  experimental  data  for  these  two  configurations. 


11.3  Description  of  Models  Used  for  Comparison 
11.3.1  Laminar  Chemistry  and  Eddy  Breakup  Models 

The  lowest  order  of  approximation  for  turbulent  reaction  rates  is  the  use  of  laminar  chemistry, 
where  the  mean  temperature  and  mean  species  mass  fraction  are  used  to  calculate  the  reaction 
rate  coefficients.  This  simplistic  treatment,  though  known  to  produce  erroneous  results,  is  still 
occasionally  used.  Because  of  the  temperature  and  species  concentration  fluctuations  in  a 
turbulent  flame,  the  reaction  rates  calculated  using  laminar  chemistry  model  are  often  too  high.  A 
quick  fix  to  this  problem  is  found  in  the  eddy  breakup  model,  where  an  upper  limit  is  set  for  the 
reaction  rate  based  on  the  turbulent  time  scale. 


11.3.2  Assumed  PDF  Model 

The  next  level  of  approximation  in  the  hierarchy  of  turbulent  combustion  models  is  the  assumed 
pdf  model.  The  assumed  pdf  method  (see,  e.g.,  [11-2])  relies  on  assumptions  of  the  shape  of  the 
pdf,  and  can  only  handle  simple  chemistry  models  where  a  single  conserved  scalar  can  be 
defined  from  which  all  other  variables  can  be  derived.  The  advantage  of  this  method  is  its 
simplicity,  while  the  disadvantages  are  its  limited  applicability  and  the  a  priori  assumption  of  the 
shape  of  the  pdf. 


11.3.3  The  Scalar  PDF  Method 

An  alternative  to  assuming  the  shape  of  the  pdf  is  to  solve  for  it  from  a  pdf  transport  equation. 
Hence  a  model  for  treating  turbulent  combustion  is  not  needed.  A  pdf  transport  equation  can  be 
derived  for  the  scalars  such  as  the  species  mass  fractions  and  the  enthalpy,  and  from  the  transport 
equation  the  pdf  distribution  can  be  solved  for;  this  forms  the  basis  for  the  scalar  pdf  method 
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[11-3,  11-4].  The  scalar  pdf  method  offers  the  advantage  of  treating  the  turbulent  reaction  rate 
without  modeling.  It  relies  on  a  finite  difference  or  finite  volume  flow  formulation  to  provide  the 
mean  velocity  and  turbulence  time  scale,  and  turbulent  convection  of  scalars  has  to  be  modeled. 
The  other  advantage  of  the  scalar  pdf  is  that  it  can  be  used  in  conjunction  with  existing,  well 
established  numerical  methods  for  the  flowfield  calculations.  However,  the  inherited  deficiencies 
of  the  current  Reynolds-averaged  methods  in  calculating  turbulent  transport  are  the  drawbacks. 
While  the  interaction  of  turbulence  on  chemistry  is  accounted  for  in  the  scalar  pdf  method,  the 
reverse  interactions  (see  section  2.1.4)  are  not  fully  accounted  for. 

The  joint  pdf  method  overcomes  the  problems  with  calculating  turbulent  transport  and  more 
comprehensively  accounts  for  turbulence/chemistry  interactions. 


11.4  Description  of  the  Solution  Methods  Used 
11.4.1  The  Scalar  PDF  Method 

A  node-based  Monte  Carlo  method  is  used  in  the  present  study  to  solve  the  scalar  pdf  equation 
[11-3,  11-4].  The  mean  flow  variables,  including  the  mean  velocity,  pressure,  and  turbulent 
kinetic  energy,  are  provided  by  a  pressure  based  variable-density  finite-volume  flow  solver 
(which  will  be  further  described  in  the  next  section).  The  species  mass  fraction,  temperature,  and 
density  are  obtained  from  the  Monte  Carlo  pdf  solver.  To  couple  the  finite  volume  solver  and  the 
Monte  Carlo  pdf  solver,  the  mean  velocity  and  turbulent  time  scale  are  passed  from  the  finite 
volume  solver  to  the  Monte  Carlo  solver  and  the  mean  density  is  passed  from  the  Monte  Carlo 
solver  to  the  finite  volume  code. 


11.4.2  Finite-Volume  Method 

The  finite  volume  method  used  for  the  scalar  pdf  calculation  as  well  as  conventional  model 
calculations  is  the  well  known  pressure-based  low-speed  variable-density  Navier-Stokes  method 
(SIMPLER)  [11-5],  with  a  standard  k-e  model.  The  solver  employs  a  staggered  grid,  with  the 
scalar  quantities  located  at  the  centers  of  the  cells  and  the  velocities  at  the  cell  boundaries.  A 
power-law  differencing  scheme  is  used  to  discretize  the  governing  equations.  Various  chemistry 
models  and  spray  and  evaporation  models  are  available  in  the  code. 
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An  assumed  pdf  model  is  used  in  the  reaction  calculations.  The  transport  equations  of  the 
mixture  fraction,  «j)>,  and  mixture  fraction  variance,  <(|)'2>,  are  solved.  A  beta-function  pdf  is 
determined  based  on  «j»  and  «|)'2>,  from  which  the  temperature  and  density  are  determined. 

When  using  multistep  chemistry  models  in  the  finite  volume  method,  the  above  assumed  pdf 
method  is  not  applicable,  and  the  eddy  breakup  model  is  used  to  simulate  the  effect  of  turbulence 
on  the  reaction  rate.  In  this  method,  the  reaction  rate  is  calculated  using  laminar  chemistry,  but  an 
upper  limit  for  the  reaction  rate  is  determined  based  on  the  turbulence  time  scale. 


11,4.2  The  Joint  PDF  Method 

The  joint  pdf  method  used  is  the  same  as  that  described  in  Chapter  10.  The  new  pressure 
algorithm  PA  II  along  with  the  SFM II  and  SLM II  models  for  frequency  and  velocity  are  used. 


11.5  Thermochemistry 

Hydrogen  and  methane  flames  are  studied  in  the  present  work.  A  fast  chemistry  model  is  used 
for  the  hydrogen  flame  calculations  because  the  time  scale  for  hydrogen-air  reaction  is  very 
small  compared  to  the  turbulent  time  scale. 

For  methane  flame  calculations,  a  general  2-step  chemistry  by  Westbrook  and  Dryer  [11-6]  for 
saturated  hydrocarbon  fuels  is  used.  The  two  steps  are: 

C„Hm  +  ^+>2  =  ”CO+fH2° 

co+fo2  =  co2 

with  the  reaction  rate  constants  for  methane  given  by 
Ai=1.5xl07,  Ni=0,  Ei=1.25xl05 

A2=1.0xl012-35,  N2=0,  E2=1.67x105  • 

For  details  on  the  chemistry  models,  the  readers  are  referred  to  Ref.  1 1-6. 
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For  both  the  fast  chemistry  and  the  2-step  chemistry  models,  lookup  tables  were  created  to 
reduce  the  CPU  requirements  of  the  calculations.  In  the  case  of  the  fast  chemistry,  a 
one-dimensional  table  is  created.  For  the  two-step  chemistry,  a  three-dimensional  table  is 
generated  for  a  given  specific  time  increment.  In  all  the  table  generation  processes,  the  NASA 
CEC  thermal  data  were  used  to  calculate  the  variable  specific  heat  and  the  temperature. 


11.6  Results  and  Discussion 


11.6.1  Swirling  Hydrogen  Jet  Diffusion  Flame 

For  the  purpose  of  quantitatively  assessing  the  various  models  used  in  combustor  designs,  the 
first  case  considered  in  the  present  study  is  a  hydrogen  swirl  combustor  [11-7].  The  schematic 
and  the  details  of  the  combustor  and  experiments  are  presented  in  Chapter  10.  The  inlet  average 
mean  velocities  are  100  m/s  for  the  central  fuel  jet,  20  m/s  for  the  annular  swirling  air  jet,  and  4 
m/s  for  external  air.  The  swirl  angle  is  30  degrees,  corresponding  to  a  swirl  number  of  0.38. 

The  hydrogen  swirling  jet  experiment  is  one  of  the  very  few  cases  where  there  is  no  ambiguity 
about  the  inflow  conditions  for  the  velocity  field,  which  makes  it  ideal  for  model  validation  and 
evaluation.  Velocities  and  turbulent  kinetic  energy  at  1.5  mm,  or  0.17  fuel  jet  diameters,  down 
stream  of  the  jet  exit  are  provided  by  the  experiment  and  are  used  in  the  present  study  as  inflow 
boundary  conditions.  The  inflow  profile  for  the  fuel  mass  fraction  is  assumed  to  be  uniform. 

The  computation  assumes  the  flow  field  to  be  axisymmetric.  For  the  scalar  pdf,  assumed  pdf,  and 
laminar  chemistry  computations,  the  computational  domain  extends  280  mm  downstream  from 
the  jet  exit;  the  lateral  boundary  of  the  domain  is  75  mm  from  the  centerline.  A  70x76  grid  is 
used  in  the  computation.  To  enhance  the  resolution  in  regions  of  interest,  the  grid  is  slightly 
stretched.  For  the  Monte  Carlo  simulation  of  the  scalar  pdf,  50  notional  particles  are  used  for 
each  computational  cell.  A  time  averaging  is  performed  to  reduce  the  statistical  error  [11-8].  The 
time  averaging  is  restricted  to  about  the  last  300  time  steps. 

In  the  results  presented  below,  the  velocity  is  nondimensionalized  using  the  centerline  nozzle  exit 
velocity  <U>OC=130.3  m/s.  The  axial  distances  and  radius  are  normalized  by  the  nozzle  diameter, 
D=9.45mm  and  nozzle  radius  R=D/2,  respectively.  The  temperature  is  normalized  by  the 
stoichiometric  temperature,  Tst=2377K. 
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Figures  11-1  through  11-3  show  the  predicted  velocity  and  turbulent  kinetic  energy  distributions 
compared  with  experimental  data.  The  symbols  represent  the  experimental  data  and  curves  are 
numerical  predictions  using  joint  velocity- scalar  pdf,  scalar  pdf,  assumed  pdf,  and  laminar 
chemistry,  as  indicated  by  the  figure  legends.  The  joint  pdf  method  predicted  the  <U>  velocity 
better  than  the  other  three  methods.  It  also  could  be  argued  that  the  joint  velocity- scalar  pdf  does 
a  better  job  for  turbulent  kinetic  energy  at  some  locations.  The  solutions  from  the  scalar  pdf,  the 
assumed  pdf,  and  the  laminar  chemistry  for  the  velocity  distribution  seem  to  overlap  one  another. 
The  joint  velocity-scalar  pdf  also  provides  information  on  other  statistical  quantities  and  higher 
moments,  that  are  discussed  in  Chapter  10. 

Of  more  importance  to  combustion  analysis  than  the  velocity  distribution  is  the  temperature  rise 
due  to  combustion.  Figure  11-4  shows  the  numerical  solution  for  temperature  compared  to  the 
experimental  data,  and  Figure  1 1-5  shows  the  joint  velocity- scalar  pdf  and  scalar  pdf  predictions 
of  the  temperature  fluctuation  compared  with  experimental  data.  One  immediately  notices  that 
the  solution  obtained  using  a  laminar  chemistry  is  unacceptable:  The  peak  temperature  is  grossly 
overpredicted,  and  the  flame  is  too  narrow.  The  assumed  pdf  method,  because  it  took  into 
consideration  the  effects  of  turbulent  fluctuations,  provided  a  much  better  solution  for  the 
temperature  field,  although  it  still  underpredicted  the  flame  edge  temperature  at  x/D=7.94  and 
overpredicted  the  temperature  at  x/D=26.5.  Both  the  scalar  pdf  and  the  joint  pdf  methods 
predicted  the  temperature  distribution  fairly  accurately,  with  the  joint  pdf  method  matching  the 
data  slightly  better,  demonstrating  the  superiority  of  the  pdf  methods. 

Figure  11-6  shows  the  temperature  variance  predicted  by  the  scalar  and  joint  pdf  methods 
compared  with  experimental  data.  The  measurement  of  temperature  variance  in  the  experiment 
depended  on  the  existence  of  nitrogen,  and  therefore  missed  the  inner  peak  corresponding  to  the 
large  temperature  gradient  on  the  left  side  of  the  flame  for  x/D=1.06,  5.29,  and  7.94.  Although 
some  differences  exist  between  the  predictions,  both  the  scalar  pdf  and  the  joint  pdf  methods  are 
able  to  predict  the  temperature  fluctuations  with  reasonable  accuracy. 


11.6.2  Step-Swirl  Methane  Combustor 

The  second  case  considered  in  this  study  is  a  step-swirl  methane  combustor  for  which 
experimental  data  were  obtained  by  Durbin  et  al.  [11-9].  The  combustor  details  are  presented  in 
Chapter  10.  The  case  considered  has  an  outer  swirl  angle  of  30  degrees  and  no  inner  swirl.  The 
average  velocity  of  the  inner  air  jet  is  14.4  m/s;  the  average  velocity  of  the  outer  air  jet  is  8.6  m/s; 
and  the  average  fuel  jet  velocity  is  2.5  m/s. 
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Figure  1 1-3.  Computed  turbulent  kinetic  energy  Figure  1 1-4.  Solutions  for  the  mean 
profiles  compared  with  experimental  data  for  temperature  compared  with  experimental  data 
the  swirl  hydrogen  flame..  for  the  swirl  hydrogen  flame.. 


The  numerical  setup  of  the  step-swirl  combustor  calculation  is  similar  to  that  of  the  swirling 
hydrogen  flame  calculation,  with  the  same  number  of  grid  points  and  notional  particles.  The 
computational  domain  extends  75  mm  from  the  centerline  and  130  mm  downstream  from  the  jet 
exit. 

The  flow  field  in  the  step-swirl  combustor  consists  of  two  recirculation  zones,  one  caused  by  the 
step,  and  a  second  one  caused  by  the  outer  swirling  jet.  The  flow  is  fully  elliptic,  and  the  ability 
to  resolve  the  recirculation  is  crucial  in  the  computation  of  this  case. 

Unfortunately,  the  experimental  data  for  inflow  conditions  is  incomplete  for  this  case:  velocity 
profiles  are  available  at  3  mm  downstream  of  the  inlet  only  from  r/R=0  to  r/R=1.3,  and 
temperatures  or  fuel  mass  fraction  distributions  are  not  available.  Because  the  exit  velocity 
profile  of  the  swirler  is  far  from  uniform,  the  information  on  the  inflow  condition  is  crucial  to 
the  correct  simulation  of  swirling  jets.  In  the  present  study,  the  inflow  velocity  profiles  had  to  be 
judiciously  extended  beyond  r/R=1.3  based  on  reported  experimental  mass  flow  rates  from 
different  streams.  With  guidance  from  the  conditioned  (conditioned  upon  different  streams) 
velocity  measurements,  a  sharp  peaked  mixture  fraction  distribution  is  used  instead  of  a  plug 
flow  profile. 

In  the  results  presented  below,  the  velocity  is  normalized  using  the  centerline  nozzle  exit  velocity 
<U>0C=21.6  m/s.  The  axial  distances  and  radius  are  normalized  by  the  nozzle  diameter, 
D=40mm  and  nozzle  radius  R=D/2,  respectively.  The  temperature  is  normalized  by  the 
stoichiometric  temperature,  Tst=2272K. 

Figures  11-6  through  11-8  show  the  comparison  between  the  predicted  velocity  profiles  and  the 
experimental  data,  and  Figure  11-9  shows  the  comparison  for  turbulent  kinetic  energy.  The 
velocity  results  show  that  the  finite  volume  method  (for  both  scalar  pdf  and  eddy  breakup  model) 
predicted  a  higher  spreading  rate  for  the  outer  jet,  or  a  larger  separation  bubble,  than  the  joint  pdf 
method.  In  general,  both  the  velocity  components  and  the  turbulent  kinetic  energy  are  well 
predicted  by  the  three  different  methods. 

The  comparison  between  the  numerical  solutions  for  temperature  and  experimental  data  is  shown 
in  Figure  11-10.  The  results  for  x/D=0.5  show  that  the  joint  pdf  method  predicted  the  convection 
of  scalars  more  accurately,  while  the  other  two  methods  predicted  too  much  spreading.  The 
advantage  of  the  pdf  methods  is  highlighted  again  in  the  temperature  predictions.  Although  in  the 
eddy  breakup  model  calculation,  the  reaction  rate  is  modified  to  take  into  account  of  turbulent 
effects,  because  the  model  does  not  account  for  temperature  and  species  fluctuations,  the  solution 
show  unrealistically  high  peak  temperature  at  x/D=0.5,  similar  to  the  laminar  chemistry  solution 
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in  the  hydrogen  case.  The  peaks  exist  irrespective  of  the  changes  made  in  the  inflow  conditions. 
On  the  other  hand,  the  pdf  models  were  able  to  correctly  simulate  the  physics  in  the  turbulent 
flame,  including  fluctuating  temperature  and  its  effects  on  reaction  rates,  the  mean  temperature  is 
better  predicted. 

The  temperature  distributions  predicted  by  the  three  different  methods  at  downstream  locations 
show  some  differences.  Although  the  predicted  temperature  profile  at  x/D=1.5  gives  the 
appearance  that  the  eddy  breakup  model  produces  an  accurate  solution,  an  overall  evaluation 
shows  that  the  eddy  breakup  model  overpredicts  the  temperature  near  the  jet  exit  while 
underpredicts  it  at  downstream  locations  (x/D=2.0),  and  the  agreement  at  x/D=1.5  is  quite 
coincidental. 

Because  of  the  uncertainties  in  inflow  conditions  and  the  sensitivity  of  the  solutions  to  these 
inflow  condition,  more  detailed  comparisons  cannot  be  made.  However,  based  on  the  results 
obtained,  one  can  draw  the  following  conclusions: 

(1)  The  eddy  break  up  model  can  produce  spurious  peaks  in  the  mean  temperature  in  the 
reaction  zone.  The  error  in  temperature  could  be  as  high  as  T/Tst=0.25  (570K)  or  higher  (see 
values  at  r/R=0.7  and  r/R=1.5~1.7,  x/D=0.5  in  Figure  11-10).  Such  large  errors  in  the  high 
temperature  regions  have  serous  implications  on  the  predictions  of  combustor  performance  and 
emissions  predictions. 

(2)  Although  the  scalar  pdf  method  accounts  for  turbulence  fluctuation  effects,  inaccurate 
prediction  of  the  mean  and  turbulent  transport  of  scalars  can  seriously  undermine  its  accuracy,  as 
can  be  observed  from  results  at  r/R>1.5  at  x/D=0.5  in  Figure  1 1-10. 

It  is  well  known  from  experimental  observations  and  detailed  chemical  kinetic  calculations  that  a 
difference  of  50K  in  mean  temperature  in  the  vicinity  of  the  reaction  zone  can  cause  large 
differences  in  pollutant  formation,  especially  for  NOx.  The  present  results  underscore  the 
importance  of  the  pdf  methods  in  accurately  predicting  combustion  efficiency  and  pollutant 
emissions  such  as  NOx,  CO  and  unbumed  fuel.  Similarly,  the  ability  of  the  joint  pdf  method  to 
accurately  calculate  turbulent  transport  is  also  very  important  for  predicting  the  above  mentioned 
combustor  parameters. 

All  the  computations  presented  in  this  work  are  done  on  IBM  RS6000  work  stations.  Because 
there  is  no  uniform  criteria  for  convergence,  a  comparison  on  exact  CPU  hours  is  not  possible.  In 
general,  the  scalar  pdf  method  is  estimated  to  be  4  to  5  times  more  time  consuming  and  the  joint 
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Figure  11-6.  Computed  axial  velocity  profiles  Figure  11-7.  Computed  radial  velocity  profiles 
compared  with  experimental  data  for  the  compared  with  experimental  data  for  the 
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Figure  1 1-10.  Computed  mean  temperature  profiles  compared  with  experimental  data  for  the 

step-swirl  methane  combustor. 
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pdf  method  5  to  6  times  more  time  consuming  than  the  finite  volume  schemes.  Considering  that 
the  pdf  algorithms  are  new  and  no  significant  efforts  are  made  to  speed  up  the  computation,  we 
believe  that  considerable  speedup  can  be  achieved  through  better  choices  of  numerical 
parameters  and  improvements  in  solution  algorithms.  Furthermore,  the  Monte  Carlo  method  is 
very  amenable  to  parallel  processing,  and  considerable  speed  up  can  be  achieved  through  parallel 
computing.  It  is  estimated  that  a  realistic  combustor  design  calculation  can  be  achieved  with  the 
joint  pdf  method  with  a  1-day  turnaround  time. 


11.7  Computations  for  a  Low  Emissions  Lean  Premixed  Combustor 

The  unique  advantages  offered  by  the  pdf  method  are  demonstrated  through  an  application  to  the 
design  of  one  of  Allison’s  advanced  low  emission  combustors.  The  combustor  is  based  on  a 
premixer  module  concept  in  which  fuel  (natural  gas)  and  air  are  mixed  under  high  swirl  in  a 
premixing  module.  Combustion  occurs  in  a  lean  premixed  mode  so  that  the  emissions  are 
significantly  lower  than  in  conventional  diffusion  flame  combustors.  A  schematic  of  the 
premixer  module  is  shown  in  Figure  11-11.  For  the  purpose  of  demonstration,  the  scalar  pdf 
method  is  used  for  the  computations  of  the  mixing  and  reacting  flow  within  the  module  as  well 
as  in  the  combustor  downstream.  (More  work  is  needed  with  the  joint  pdf  method  for  application 
to  complex  geometries  with  body  fitted  grids  as  explained  in  the  next  chapter.  However,  the 
advantages  seen  for  the  scalar  pdf  method  will  only  by  enhanced  with  the  joint  pdf  method.  The 
two-step  finite  rate  methane  chemistry  (section  11.5)  is  used  for  the  computations.  Figure  11-12 
shows  the  contour  plots  of  the  computed  mean  axial  velocity  field,  along  with  those  of  mean 
temperature  and  fuel  mass  fractions.  Measurements  for  mean  velocities,  mixture  fraction  and 
mean  temperature  have  been  made  and  data  are  available  for  comparison.  The  computations  are 
in  good  agreement  with  data  both  qualitatively  and  quantitatively. 

The  crucial  advantage  provided  by  the  pdf  method  is  the  ability  to  more  realistically  model 
premixed  combustion  with  the  Monte  Carlo  algorithm.  The  joint  velocity-scalar  pdf  method  has 
been  used  previously  in  detailed  modeling  studies  of  premixed  flames  [11-10  through  11-12]  and 
its  ability  to  calculate  turbulent  flames  speeds  and  other  fundamental  parameters  has  been 
demonstrated. 

For  the  combustor  under  consideration,  the  results  in  Figure  11-12  show  that  the  flame  is 
stabilized  in  front  of  the  large  recirculation  zone  as  expected  since  the  velocities  are  low  in  that 
region.  The  extent  of  the  recirculation  zone  and  the  location  of  the  flame  are  consistent  with  the 
experimental  data.  Also,  as  expected  and  borne  out  by  experiments,  the  discrete  fuel  jets  are  well 
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mixed  within  the  premixer  module  and  the  mixture  fraction  distribution  at  the  exit  of  the 
premixer  module  and  especially  in  the  combustor  is  uniform.  Since  the  mixture  fraction  within 
most  of  the  module  is  lean  and  within  combustible  limits,  conventional  finite-volume  based 
combustion  models  predict  the  flame  to  spread  well  into  the  premixer  model  (even  with  using 
finite-rate  chemistry)  thereby  predicting  "flashback"  as  well  as  erroneous  temperatures  and 
emissions.  Accurate  prediction  and  avoidance  of  flashback  is  one  of  the  key  issues  in  premixed 
combustor  design.  However,  with  the  pdf  method  used,  even  if  the  mixture  fraction  is  within 
combustible  limits,  the  individual  computational  particles  do  not  burn  until  they  come  in  contact 
with  other  burnt  particles.  Hence,  the  interaction  of  convective  and  turbulent  transport  with 
flame-speeds  is  captured  resulting  in  the  correct  prediction  of  the  flame. 


Figure  11-11.  Schematic  of  the  RSPN  low  emissions  combustor  module. 
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Premixer  Module  Combustor 
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Figure  11-12.  Contour  plots  of  mean  mixture  fraction,  temperature  and  axial  velocity 
computed  using  the  PDF  method  for  the  RSPN  low  emissons  combustor.  PDF  methods 
offer  unique  capability  for  premixed  combustor  design.  (Scalar  PDF  is  used  here.) 


11.8  Summarizing  Remarks 


Five  different  combustion  models  were  evaluated  in  the  present  study.  The  results  show  that  the 
use  of  laminar  chemistry  produces  unacceptable  results,  and  the  eddy  breakup  model  also 
produces  inaccurate  temperature  predictions.  Assumed  pdf  is  a  major  improvement  over  the 
other  conventional  combustion  models  and  provides  reasonably  good  aerothermal  predictions, 
provided  that  simple  one  step  chemistry  models  are  acceptable.  The  present  work  again 
demonstrated  that  the  scalar  pdf  and  joint  velocity-scalar  pdf  methods  are  superior  to  the 
conventional  combustion  models  in  handling  turbulent  combustion  problems.  The  present  results 
underscore  the  importance  of  the  pdf  methods  in  accurately  predicting  combustion  efficiency  and 
pollutant  emissions  such  as  NOx,  CO  and  unbumed  fuels.  This  is  demonstrated  by  an  application 
to  a  "real"  combustor  hardware.  The  additional  ability  of  the  joint  pdf  method  over  scalar  pdf 
method  to  accurately  calculate  turbulent  transport  and  other  turbulence/chemistry  interactions  is 
also  important  for  predicting  combustor  performance  and  emissions. 
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12.  CONCLUDING  REMARKS 


Significant  progress  has  been  made  in  the  development  of  an  advanced  and  more  accurate 
combustor  design  tool  based  on  the  joint  probability  density  function  (pdf)  method.  The  pdf 
method  was  chosen  based  on  its  inherent  and  demonstrated  advantages  over  conventional  other 
advanced  methods  for  treating  turbulent  reacting  flows. 

Several  new  submodels  and  algorithms  within  the  pdf  method  were  developed  and  validated  to 
enable  the  application  of  the  pdf  method  to  complex  flows  of  practical  interest,  such  as  in  gas 
turbine  combustors.  The  method  and  the  models  have  been  validated  against  a  variety  of  flows 
consisting  of  the  essential  features  of  gas  turbine  combustor  flows.  Both  available  experimental 
data  and  experiments  specifically  designed  to  fill  the  voids  in  data  in  the  literature  were  used  to 
provide  benchmark  quality  data  for  the  validation.  The  joint  pdf  results,  were  uniformly  in  good 
agreement  with  data  for  all  the  flows  considered.  The  joint  pdf  results  were  compared  with 
results  from  conventional  methods  and  assessed  against  experiment  data.  The  joint  pdf  results 
were  clearly  more  accurate  than  those  from  conventional  methods.  A  baseline  combustor  design 
system  based  on  the  joint  pdf  method  has  been  developed. 

However,  significant  development  effort  lies  ahead.  There  are  several  areas  for  further 
improvement  of  the  baseline  design  system  so  that  the  full  potential  of  the  pdf  design  system  is 
realized.  The  development  of  the  method  under  the  current  program  focused  on  gaseous  fuels. 
The  current  design  system  can  be  applied  to  gaseous  fuel  applications  such  as  in  certain  low 
emissions  combustor  concepts.  There  are  no  theoretical  limitations  for  incorporating  liquid  fuel 
sprays  in  the  pdf  method.  In  fact,  the  advanced  spray  transport  models  are  also 
Lagrangian-based  and  use  the  Monte  Carlo  technique.  Hence  the  pdf  method  and  the  spray 
transport  models  are  naturally  compatible.  A  preliminary  demonstration  of  the  pdf  method 
incorporating  sprays  (although  nonevaporating  and  nonreacting)  in  the  boundary-layer  algorithm 
has  been  performed  by  Anand  [12-1].  The  method  needs  to  be  extended  for  evaporating  and 
reacting  sprays  and  incorporated  in  the  elliptic  flow  algorithm.  Allison  is  a  leader  in  the 
development  of  advanced  evaporation  and  spray  models  [12-2,  12-3].  These  models  will  be 
incorporated  in  the  joint  pdf  design  system. 

The  other  areas  of  further  development  are  the  treatment  of  internal  blockages,  development  of 
the  pdf  method  using  body-fitted  grids,  accurate  and  efficient  implementation  of  the  detailed 
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combustion  chemistry,  enhancement  of  computational  efficiency  through  advanced  parallel 
processing  strategies,  and  extension  to  complex  3-D  flows. 

Treatment  of  internal  blockages  and  the  use  of  body-fitted  grids  are  needed  to  model  complex 
combustor  geometries  and  the  flows  in  the  combustor  diffuser  or  casing.  The  main  challenge 
here  is  to  develop  the  tracking  procedure  for  irregular  grids  in  the  Monte  Carlo  solution 
algorithm  and  to  develop  a  pressure  algorithm  for  such  grids. 

To  fully  exploit  the  advantages  offered  by  the  pdf  method  to  treat  complex  chemistry,  a  method 
to  include  detailed  combustion  (which  typically  includes  hundreds  of  reactions  and  tens  of 
species)  in  a  computationally  efficient  way  is  needed.  Allison  is  collaborating  with  Professor 
Pope  of  Cornell  University  to  further  develop  and  validate  such  a  method  which  is  in  part  based 
on  the  Manifold  Method  [12-4,  12-5].  The  method  has  shown  considerable  promise  and  has 
demonstrated  that  accuracy  matching  that  of  the  detailed  chemistry  can  be  obtained  with 
considerably  less  computational  effort. 

As  more  complex  flows  are  calculated  with  the  pdf  method,  it  is  essential  to  utilize  parallel 
processing  strategies  with  the  goal  of  achieving  computational  turnaround  times  of  a  day  to  allow 
the  combustor  designer  to  effectively  analyze  and  iterate  on  combustor  designs.  The  current 
baseline  design  system  already  incorporates  basic  parallel  processing  strategies  and  uses  particle 
partitioning.  Other  strategies  such  as  domain  decomposition  and  other  advanced  techniques  need 
to  be  explored. 


Computations  in  complex  3-D  geometries  will  also  need  to  exploit  multiblock  and  multigrid 
solution  strategies  to  reduce  computational  effort  and  to  have  geometry  flexibility  without 
sacrificing  computational  accuracy.  A  version  of  the  multigrid  capability  has  already  been 
implemented  and  tested  in  the  baseline  design  system. 

In  conclusion,  the  current  program  has  been  instrumental  in  the  significant  development  and 
validation  of  the  joint  pdf  method  as  a  gas  turbine  combustor  design  tool.  The  demonstrated 
level  of  accuracy  that  can  be  achieved  with  the  pdf  method  in  calculating  complex  chemistry, 
turbulent  transport,  and  turbulence/chemistry  interactions  is  key  to  the  accurate  prediction  of 
critical  combustion  performance  characteristics  such  as  efficiency,  emissions,  exit  and  wall 
temperatures,  lean  blowout,  stability,  flashback,  autoignition,  etc.,  in  advanced  combustor 
concepts.  Several  planned  future  development  activities  will  lead  to  the  realization  of  the  full 
potential  of  the  joint  pdf  design  system. 
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